CLUTCH SIMULATION FINAL REPORT TO FORD MOTOR CORPORATION by C. Pierre, Assistant Professor J.R. Barber, Professor P. Papalambros, Associate Professor P. Cha, Graduate Student April 1986 #UM-MEAM-86-38

CLUTCH SIMULATION FINAL REPORT TO FORD MOTOR CORPORATION by C. Pierre P. Cha Assistant Professor Graduate Student J.R. Barber Professor P. Papalambros Associate Professor UM Contract No. 388792 April, 1986

TABLE OF CONTENTS (1) Numerical Simnulation of Clutch Dynamics a. Introduction b. Current Status c. System Model d. Numerical Solution e. Programming Logic f. Program Usage g. Future Goals (2) Supplementary Tables and Diagrams a. Parameter Description b. State Variable Listing c. System Control Diagram d. Programming Hierarchy (3) Programming Logic of Clutch.LKut a. Primary Level b. Secondary Level: Interactive Options c. Differential Equation Subroutines Called by KUTTA for Appropriate State d. Forcing Function Called by KUTTA to Provide Input Data e. Data Base Subroutines Called by State Subroutines (4) Program Flow Chart (5) Graphical State Representation (6) Translational Model

a. Model Description and Parameters b. Initial Clutch Positions and Physical Clutch Parameters c. System States d. Translational Equations of Motion Summary (7) Rotational Model a. Model Description and Parameters b. Rotational Equations of Motion (8) Analysis of Coulomb Friction (9) Results a. Case 1: Clutch Engagement b. Case 2: Engine Torque Transmission (10) Program Listing

1. NUMERICAL SNMULATION OF CLUTCH DYNAMICS l.a Introduction The primary purpose of this report is to provide Ford with a clear understanding of the current status of the main clutch simulation program. Material is presented which documents the system model, numerical solution method, programming logic, and program usage. l.b CuTrent Status Major accomplishments since the last update are as follows: A new numerical integration scheme has been implemented to replace the IMSL subroutine DGEAR. For the major drawback of DGEAR was that it could not jump accurately from one set of equations to another, i.e., from state to state. The program now uses a fourth order Runge-Kutta integration scheme. Very accurate numerical results have been obtained. In addition, the system now passes from state to state while undergoing stiction/friction phenomenon. * Obtained numerical solutions have been cheched against known analytical results, and excellent agreement has been observed. Numerical results are presented in Section 9. * Previous errors associated with incorrect axes scaling in subroutine GRAPH have also been eliminated.

l.ce System Model The system model is described on several levels which are covered thoroughly in Sections 2 through 7. These include a System Control Diagram, and the state equation derivations with corresponding component free body diagrams. These are fully up dated and fairly self explanatory. At this point, some comments can be made concerning current limitations of the overall system model. By reviewing the System Control Diagram it can be seen that the clutch is only one of several important components. Thus far, other components have been treated with great simplicity. For example, the engine is modeled as a constant torque/power output device. Obviously, engine output is a complicated function of many variables. In addition, the load or drive system and vehicle inertia are made up of a few lumped components. The pedal linkage has not been modeled at all. Instead, a time dependent forcing function (Finjou) is applied directly to the release bearing. Finally, human interaction has not been accounted for in the code. Many of these deficiencies have been considered in the present code in such a way that when data becomes available, it may be inserted into the appropriate subroutines immediately. However, it may be proven necessary to increase modeling detail as well, by using a finer distribution of lumped springs, masses, and dampers. Present data used in the code refelects closely the actual clutch. An exception to this is the current linear approximation of highly non-linear characteristics. Also, it should be noted that the translational damping in the clutch model is now of viscous type. This should be changed to coulomb type. Another possibly significant effect has been temporarily neglected by not including momentum transfer and impact. This may be included between translational components such as the pressure plate and clutch assembly, or rotational components such as the outer clutch hub and inner clutch hub (transmission input shaft). Another area requiring enhancement is the stiffness and kinematics of the bellville spring, clutch cover, and drive straps. These are now lumped into one spring and mechanical lever. While this may be acceptable from a quasi-static point of view (i.e. stiffnesses may be added), it may not hold for a dynamic model. Each component should be modeled independently Also, because of initial flex in the clutch cover during disengagement, the bellville spring

actually pivots about two different points; first it pivots at the pressure plate, then it pivots at the clutch cover. Therfore, there is a need for defining an additional translational state and logic to determine the transition point within KUTTA, the subroutine which determines the correct system state. Lastly, it remains to apply dependent friction coefficient logic and data. This includes temperature effects, velocity effects, and static/slipping consideration. There can be little question that the friction coefficient treatment is very important to obtaining an accurate model. A current error in the code is the inability of the torsional friction damper in the clutch hub interface to resist motion in either direction. The same is true at the rear wheel. That is, a change in sign on the friction term is necessary in the equation of motion which describes the appropriate degrees of freedom where coulomb friction is concerned. Otherwise, energy ill be incorrectly input to the system. This problem has been addressed in the main clutch pad/lpressure plate/flywheel interface (see flow chart KUTTA). Since there are so many friction dependent equations of motion, it might make sense to create a seperate, generalized friction/stiction subroutine which could be routinely called. This would ease the modeling expansion process greatly.

1.d Numerical Solution The numerical solution method used in the earlier version of the program was an integration subroutine called DGEAR. This subroutine solved a simultaneous system of first order differential equations. It had the option of either using the Adams solution method or the Gear solution method (for stiff system). Essentially, DGEAR was called at the incremental time steps, through a DO-loop in the main program. Each time DGEAR was called, it in turn called subroutine FCN. FCN was responsible for determining the correct system state, whereupon that state was then called. Once this had occured, the state variable first order derivatives were calculated and passed back through FCN to DGEAR. At this point DGEAR was able to carry out the integration process to determine the appropriate state variable values which could then be used as input for the next time step. The major drawback of this method is that DGEAR cannot jump from one set of equations to another. This is due to the fact that DGEAR does not need to call FCN at every time step. DGEAR is a very efficient method because of its ability to extrapolate future values. It then corrects itself upon the next call to FCN. It, however, causes a transition between state equations to be delayed since FCN is only called when DGEAR deems it necessary. This problem was found to be very significant, since it was observed that the use of DGEAR led to erroneous results. Thus, the program had to be modified by using a different solution method, namely the fourth order Runge-Kutta time integration scheme. The principal advantages of the Runge-Kutta method lie in its self starting features and its ease of programming. Here, it also has the added advantage in that it enables a smooth transition between state equations to occur since subroutine KUTTA (the new subroutine which determines the correct system state) is called at every time step. Thus, if the time step is small enough, transition between states can be effectively and quite accurately simulated. A fourth order Runge-Kutta method employs a recurrence formula of the form (k, +2k+2k3+ k,4) 6 to calculate successive values of the dependent variable y of the differential equation -d= f(y t)

The k,'s are as follows, k ALS f(ti, Yi) = A(ti + t,y y) k3 = Atf(ti + Y, yi + k4 = f(tit + A, Yi +'3) In this fourth order Runge-Kutta method, the per step error is of the order (At)5. Subroutine KUTTA has the same function as FCN in the earlier version. It determines the correct system state, then calls that corresponding state. Once this has occured, the state variable first order derivatives are computed and passed back to KUTTA, at which point KUTTA may perform the Runge-Kutta time integration process to determine the appropriate state variable values which can then be used as input for the next time step. 7~

L.e Programming Logic The code is of modular design for ease of revision and understanding. The section contains a programming hierarchy chart and a brief written description of each subroutine to accompany it. Of further use are parameter and state variable listings. These sources of information should be combined for a full interpretation. Essentially, the code hierarchy works as follows: First, the main program, or primary level, is run interactively. There are several secondary subroutine options which may be selected from this point. They include DEFALT, MODIFY, and GRAPH, all of which are for pre- and post-processing purposes such as parameter selection, modification, and graphics. As described in the numerical solution section, the Main program calls KUTTA. Subroutine KUTTA's primary prupose is to select the appropriate state equations, evaluate them using fourth order Runge-Kutta time integration method, and then send information back to the main program. KUTTA must therefore contain the friction/stiction logic previously discussed in order to choose the correct state. KUTTA always calls FINPUT and TORQIN. These subroutines provide forcing functions which act on the state equations. These are the release bearing forces application and engine torque output. Generally speaking, these subroutines require input information since it is a closed-loop system. Once KUTTA has selected the corrected system state, database subroutines are called to evaluate system parameters such as spring rate and damping rates. The reason for seperate subroutines is that these parameters may be highly non-linear, and therefore depend on current state variable values. these subroutines are BELSPR, MARCEL, and ENGINE. Notice that additional subroutines for friction coefficients and perhaps new parameters associated with a refined model may be easily appended to the code. l.f Program Usage The program is designed to be used interactively. Within the interactive menus are examples which explain their usage. There are currently three basic options. There are preprocessing, running the code, and post-processing. The pre-processing includes the ability to modify parameter values or select default values. Eventually, this vwill be expanded to enable initial values of state variables to be altered. Code running options include time step size and 8

total number of time steps for the system response. The post-processor contains a graphics package which plots any chosen variable versus time. 9

l.g Future Goals Future goals have been mentioned throughout this report. They shall be listed here as follows: (1) Develop a friction coefficient subroutine for the clutch/flywheel interface to take into account static and slipping behavior, as well as temperature and velocity dependence. (2) Consider the addition of a general stiction/friction logic subroutine. (3) Correct all potential energy input situations which will not occur in reality. These include angular and translational momentum conservation, impact losses, and friction force opposition to forward motion. (4) Enforce the clutch process to be reversible. Presently, the code is intended for engagement only. Universal application to various situations might necessiate some further programming logic. Work is underway concerning this topic. (5) Implement a two stage spring and possible backlash effects into the clutch hub. (6) Replace temporary numerical values in code with parameters. In general, these are currently associated with the bellville spring linear approximation and certain displacement constants. (7) Change viscous damping in clutch model to Coulomb damping. (8) Replace the single spring of stiffness k2 with independent components for the beliville spring, clutch cover, and drive straps. Also, add an additional translational state which would correctly model the tuo-pivot action of the bellville spring. (9) Include accurate engine mapping information in subroutine TORQIN. (10) Determine representative release bearing load versus time curves for use in subroutine FORCIN. (11) Refine drive train model to include more springs, masses, dampers, etc. Consider backlash. (12) Model the hydraulic/mechanical pedal linkage. (13) Correct subroutine MODIFY to work interactively. (14) Add new options to interactive menus as described in program usage section.

Obviously, these goals must be attemped in a correct priority. This topic warrants further discussion with Ford. 11

2. SUPPLEMENTARY TABLES AND DIAGRAMS 2.a Parameter Description rn-pressure plate mass m4-clutch assembly mass I1 —pressure plate/flyAheel/crankshaft/clutch cover moment of inertia 12 —outer clutch pad assembly moment of inertia 13-transmission moment of inertia 4 —vehicle mass transformed into equivalent moment of inertia Ie —engine block moment of inertia k,-engine mount sprint rate k2- release bearing load spring rate k3-marcel (cushion) spring rate ce —engine mount viscous damping rate c2 —release bearing viscous dampring rate c3 —marcel (cushion) viscous damping rate T,-engine output torque F,,,,t —release bearing input force Torqin —engine torque which drives system inner and outer clutch THKT-temporary constant used to calculate release bearing load due to spring k2 rD-radius of driving wheel rH-mean radius to clutch hub damping surface rc —mean radius to clutch pad torque transfer surface rE —radius from engine crankshaft axis to engine block mount u/ —clutch pad friction coefficient FD-tire rolling friction force 12

T2 —clutch pad torque transfer from engine M —normal force on pressure plate side clutch pad 1\2 —normal force on flywheel side clutch pad A'3-normal force on clutch hub coulomb friction damper TR-transmission gearing ratio a —ratio of release bearing displacement to pressure plate displacement Dt-initial gap between clutch pad and flywheel plus compressed clutch width 13

2.b State Variable Listings Displacements Velocities Accelerations Q(1)=z —-release bearing QDOT(1)=Q(2)=i2 QDOT(2)=i 2 (translation) Q(3)=z4-clutch assembly QDOT(3)=Q(4)=i4 QDOT(4)=Z4 (translation) Q(5)=, —engine block QDOT(5)=Q(6)=6 ~ QDOT(6)=, e (rotation) Q(7)=8, -crankshaft/flywheel QDOT(7)=Q(8)= 1 QDOT(8)='1 (rotation) Q(9)=62 —clutch pad assembly QDOT(9)=Q(10)=62 QDOT(10)=62 (rotation) Translational Clutch Model Rotational Clutch Model.ll l1'Kt e

2.c System Control Diagram TORQUE FEEDBACK TORQUE FEEDBACK TORQUE IN I TORQUE OUT ENGINZE CLUTCH LOA RPM IN RPM OlUT THROTTLE CLUTCH PEDAL POSMON FORCE HUMAN CONTROLLER ENGINE RPM VELOCITY ACCELERATION JERK 15

2.d Programing Hierarchy DEFALT MODEFY KT PRINTA GRAPH STATE 1 tc a STATE4,I I[ STATE2 SAE JJJ-2 DATA $TATE3 TATE4C DATA JJJ = 0 STATE3C TORQIN DATA B PR MARCEL ENGINE JJJ is a slip/stick indicator within the program. 16

3. PROGRAMMING LOGIC OF CLUTCH.KUT Clutch.Kut is the primary code being developed for the modeling of the clutch system. The code incorporates a fourth order Runge-Kutta time integration scheme to handle the numerical solution. While subroutine KUTTA performs the numerical solution, other subroutines provide the logic, mathematical functions and data bases necessary to model the system. To better understand Clutch.Kut logic, a programming hierarchy chart has been included, as well as a flow chart. A brief description follows. 3.a Primary Level NHIN - The main program performs several functions. (a) Initialization of system variables and parameters by calling subroutine DEFALT. (b) Initialization of interactive command sequence by offering user several primary options. These options include output graphics (subroutine GRAPH), modification of system variables or parameters (subroutine MODIFY), execution of program (subroutine KUTTA), and termination of program. (c) Writes time, state, and system variables Q(1-10) for each time step on file 08. (d) Writes time, system variables Q(1-10), as well as Q(11) (=Q(7)-Q(9)) and Q(12) (=Q(8)Q(10)) for each time step on file 07 for post-processing usage, namely through the use of subroutine GRAPH. 3.b Secondary Level: Interactive Options (1) MODIFY - Allows user to set system parameters to default values by calling DEFALT. Instead, the user may elect to alter parameters interactively by renewing appropriate constants or entering discrete data points for non-linear parameters such as a bellville spring curve. (2) DEFALT- Sets system parameters to previously designated values. (3) GRAPH - Reads dependent variable as specified by user into in y array, while reading a corresponding x value of time. Assembles data for graphing purposes. Generates 17

x-y plots with axis labels provided by user. Software incorporates *PLOTSYS plotting routines. Requires Tektronix terminal or one which emulates a Tektronix. 18

3.c Differential Equation Subroutines Called by KUTTA for Appropriate State (1) STATE1 - The first translational state of the clutch. No contact between clutch components exists here. (2) STATE2 -The second translational state of the clutch. The pressure plate and clutch face are in contact. (3) STATE3(JJJ=O) -The third translational state of the clutch. All translational clutch components are in full contact. The system is not fully compressed however. Interface surfaces are slipping. In the program, this state is designated by JJJ=O. (4) STATE3(JJJ=2) - The third translational state of the clutch. This state is identical to STATE3(JJJ=O) in the translational sense. The immediate histoy of the interface surfaces is one of non-slip. Therefore the clutch is contacting the flywheel and the pressure plate. The interface surfaces are beginning to slip. This state is designated as JJJ=2 in the program. (5) STATE3C - The third translational state of the clutch. Pressure plate and flywheel are both in contact with clutch. Interface surfaces are not slipping, they are sticking. Inner clutch components are not fully compressed. (6) STATE4(JJJ=O) -The fourth translational state of the clutch. All translational components are fully compressed. The clutch face is slipping. (7) STATE4(JJ=2) - The fourth translational state of the clutch. All translational components are fully compressed. The immediate history of interface surfaces is one of non-slip. The interface surfaces are beginning to slip. (8) STATE4C - The fourth translational state of the clutch. All translational components are fully compressed. There is no slipping at the clutch face. 3.d Forcing Function Subroutines Called by KUTTA to Provide Input Data (1) FORCIN -Supplies time dependent clutch pedal effort. In the future, this subroutine will be dependent on other variables as well, such as vehicle acceleration. This requires further insight to human response and some optimal criteria. (2) TORQIIN - Supplies engine torque input values. These values are highly influenced by 19

engine rpm., throttle setting, and load. This subroutine will be expanded to consider these items. 3.e Data Base Subroutines Called by State Subroutines (1) BELSPR - This subroutine computes values for the release bearing stiffness rate, damping rate, and hysterisis effects as functions of system variables. These rates are designated by kr and c2 respectively. Non-linearity would require the use of equations whose coefficients are determined from MODIFY and LSTSQR. This is left for future work. (2) MARCEL - Similar to BELSPR, MARCEL employs techniques which compute values of the cushion spring rate and damping factor as functions of system variables. These rates are referred to as k3 and c3 respectively and are in general, non-linear. (3) ENGINE - Computations are made for engine mouting spring rates and damping factors.

4. PROGR.AM FLOW CHART Subroutine FCN Obtain system pa ame'. Obtain system inputs. No is Yes x2>parl _ ~~~~~Is State 1 x4>par2 Return Yes Is No -—'x2>par3 State 2 cN~ No /c\ Rewmr CheckJJJ Vmew-Q(8)-Q(10) Is No VCmew)*(vrDld)4)0 it~Yes JJtate3 Compute AA Realmm No Is Yes ABS(AA )>Fr' stuck Sip, 21'{ —-

4. PROGRAM FLOW CHART(continued) lijul~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~i CheckJJJ Vmew.Q(8)-Q(10) /~I~s \No (Vmew)*(Vrowe_-nYes es ~ ~~~~~State,4.~~~~~.. Compute AA Retumn i Stuck ABS(AA)>F / I \/~~~~ slip Sta4c 1L State4 _. 22tate~e 3~JJ v2 Retum Retun 22

5. GRAPHICAL STATE REPRESENTATION contact between contact between inner clutch slipping between slipping between state first friction second friction components first friction second friction number pad and pad and fully pad and pad and pressure plate flywheel compressed pressure plate flywheel #1 NO NO #2 YES NO NO YES #3 YES YES NO YES YES #,3C YES YES NO NO NO #r4 YES YES YES YES YES #'4C YES YES YES NO NCY Note: As far as the translational equations are concerned, states #3 and #a3C are the same state. Also, translational equations do not exist for states #4 and #4C. This is due to the fact that for these states, translational displacement of the clutch components does not exist, since the inner clutch components are fully compressed. This is due to physical constraint. As far as rotational equations are concerned, states #3 and #4 are the same; states #3C and #4C are the same. 23

6. TRANSLATIONAL MODEL 6.a Model Description and Parameters 1) Pt) ulsre (1) -Pressure plate * masS=m2 * displacement=a2x2; where a2 is a displacement coefficient. (2) -First friction pad *mass=negligible ~ displacement=irrelevant (for state #1) ~ displacement= a2 (for other states) (3) — Outer and inner hub assembly * mass=m4 * displacement=z4 * coefficient of friction at interface=0.0 (4) -Second friction pad ~ mass=negligible ~ displacement=irrelevant (for all states) 24

(A) -Diaphragm spring (represented by lever A and k2 and c2) v force input=f(x2, t), where x2 is the displacement of clutch pedal. (B) -Marcel or cushion (represented by k3, C3)

6.b Initial Clutch Positions and Physical Clutch Parameters d3 N C Initial Positions Physi.cal Parameters d.=initial distance between 13=free length of k (marcel) pressure plate and cover t,=thickness of entire clutch ds=initial distance between t2=thickness of pressure plate first friction pad and flywheel t3=ts=thickness of friction pad d4=initial distance between t4=thickness of marcel clutch hub and flywheel dB=initial distance between second friction pad Initial positions are determined by user.

6.c System States State #1! -From clutch pedal fully depressed to point where pressure plate is just about to contact the first friction pad. -Engine at idle speed. State #2 -From when the pressure plate contacts the first friction pad to the point where the second friction pad is just about to contact the flywheel.

6.c System States(continued) State #3, #-3C -From when the second friction pad contacts the flywheel, to the point where the inner clutch assembly is just about compressed (i.e. clutch pedal is just about fully released) States #4, #4C -When the inner clutch pedal assembly is fully compressed (i.e. clutch pedal is fully released)

6.d Translational Equations of Motion Summary State #1 m'a2: +t c2a"i; + k2az2 = k THKT - Fin.st where THKT=constant=9.34 State #2 m.:a2" + (C2 + c3)a2i2 + (k + k3)a.22 = k2 THKT - rin,,st - &(D, - Z) + C~, max + c3i + k3; = C3a2i2 + 3(Dt + az2) where THKT=9.34, Dt=1.8

6.d Translational Equations of Motion Summary (continued) State #-3 m2 a2. + (c2 + c3)a22 + (k + k3)a2, = k THKT - Fnput - k3(Dt - z) 4 = 1.0 Mwhere THKT=9.34, Dt=1.8 THKT=Temporary constant used to calculate release bearing load due to linear spring k2. Dt=lnitial gap between clutch pedal pad and flywheel plus compressed clutch width.

Calculations of N1 and N2 State #2, #3, #3C Free-body diagram of first friction pad k3(Dt + a2x2 - z4) C3(a2Z2 - 4) 1 = k3(Dt + a.z2 - z4) + 3(a2~ - i4) State #3, #3C Free-body diagram of the second friction pad k3(Dt + a2x2- x4) C3(a2. - 4) N1 = A72

State #4, #4C Free-body diagram of hub assembly and both friction pad N2 - I N1 N1 = N2 = (:)

7. ROTATIONAL MODEL 7.a Model Description and Parameters BB. (1) — Flywheel, cover,and pressure plate * mass moment of inertia=II * angular displacement —1 (2) -Outer hub * mass moment of inertia=12 * angular displacement=82 (3) -Inner hub and transmission * mass moment of inertia=I3 (4) -Drivetrain (from output of transmission -to rear wheel) * mass moment of inertia=14 (A) -First friction pad * mass moment of inertia=negligible * angular displacement=irrelevant

(B) -Second friction pad ~ mass moment of inertia=negligible * angular displacement=irrelevant

7.b Rotational Equations of Motion: Derivation for state #1 I, _ V 11l = Te O8e: /ee + ke Iclc + cerf@ +$ kcrc35 = Tc 35

7.b Rotational Equations of Motion: Derivation for state #2 TtII 11, = T - T1 T, = pi ATl r, N, ks3(Dt + a2x2 - x4) + C3(acx ) i4) 02 T1 =,, L. 3

et I' ke Ce 1Me, + cereOe + kere e = TI

7.b Rotational Equations of Motion: Derivation for states #3 and #4 81 T,:= p, A[l r] T2= /IN2rc N1 = N2 I.il = Te - 2t/l r 82 T "T. iz (I2 + 13 + 14)82 = T: + T = 2 jN=r,

(lb f%~~~~~~~~~~~~~~~~~~~~~~~r r~~ ~ r crb It ft W ft~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~` CD +~~~~~~~~~~~~~~~~~~~~~~c,P4 m 11 ~ ~ cl I~~~b — ~~ip

7.b Rotational Equations of Motion: Derivation for states -3C and #4C 01 (I1 + 12 + 13 + )61 = Te ~Ceog~,e8 + CrcOe + kerO, = T= 40

8. ANALYSIS OF COULONIB FRICTION The model shown in Fig.1 was selected to study coulomb friction. The idealized relationship between VI (relative velocity) and Ff (friction force) shown in Fig.2 is assumed. Thus we can state that Ff = -kNsign V,. The case of V, = 0, that is sticking deserves some comment. During the motion of the masses ml and m2, two possible events can happen, either the two masses wtill stick to each other, or the two masses will slide. If the two masses start to slide, then the detection of the correct direction of Ff is important. Table 1 summarizes the procedure to detect sliding and to find the correct direction of friction force for an equivalent translational model.

F, MI * u yN lr = i1 -x Fig 1 Fig 2 Sliding or stiction F, + F2 _ F_ SIG( F + F ) In case of sliding ml + rm ml mI +m2 ml FBD for mass m sliding ml Positive ~ I >|m > m | Negative II Table 1 42

9. RESULTS Two cases are investigated and a brief description of each follows: Case 1 shows the engagement of a clutch as it goes through all the different states. The initial conditions and system constants are listed in Table 2. The system's motion occurs sequentially from State 1 to State 4. Sticking occurs in State 3. The results obtained in States 129 and 4 have been checked against exact analytical solutions, and excellent agreement has been observed. In Case 2, a constant step torque of 77000 is added to the original engine torque of 10000 at State 4. The initial conditions of State 4 are obtained from the final conditions of Case 1, and are listed in Table 2. The step torque is applied at time step number 400, or at t =.004 second. All the system constants are the same as in Case 1. M\aotion only occurs in rotational State 4. At first, the interface surfaces are stuck. When the step torque is applied, the friction force is overcome and slipping occurs. The step torque necessary to overcome the friction force has been determined theoretically and confirmed by numerical results.

Table 2 Case 1 Case 2 Initial Conditions -10.5 0.0 22 0.0 0.0 Z4 0.0 1.0 4, 0.0 0.0 8, 0.0 0.01711 8c 0.0 -0.1987 81 0.0 10.86 81 52.36 52.76 9: 0.0 3.626 ~' 0.0 52.77 System Constants Te 10000.0 10000.0 Finpu, 0.0 0.0 ka.) 200.0 200.0 C2 0.0 0.0 k3 3000.0 3000.0 C3 0.0 0.0 k, 2000.0 2000.0 Ce 10.0 10.0 44

Case 1 ~~0.00 0.04 0.08 0.12 0.16 0.20 TIME CD "' ~.00 0.04 0.08 0.12 0.16 0.20 TIME

Case 1'b.00 0.04 0.08 0.12 0.16 0.20 TIME Co cT.R - i I I i i i i i I i C.00 0.04 0.08 0.12 0.16 0.20 TIME

9 e. (X10O) e.0-].0 - 0.20 -1.00 -0.20 0.60 1.40 0.20 0.30 H 1 H 3 4 P \ FT' ()

Case 1 coO.. QO CD.00oo 0.04 0.08 0.12 0.16 0.20 TIME co o4 o., c.oo 0.04 0.08 0.12 0.16 0.20 TIME

o 0139.00 59.00 4.00 8.00 1200 Cot I I cm CD c) U. C"

Case 2 C. CD,I.00oo 0.02 0.05 0.07 0.10 0.12 TIME (X101) o).CD D</ 0.00 0.02 0.05 0.07 0.10 0.12 TIME (X101)

Case 2 (0 co ~0.00 0.02 0.05 0.07 0.10 0.12 TIME (XlO')'I CD "b.oo 0.02 0.05 0.07 0.10 0.12 TIME (X101)

Case 2 co c, C" %. Do 0.02 0.05 0.07 o.lo 0.12 TIME (X101) cd co TIME (X10') TIME (X101)

Case 2 co O C4 x (DI'T0.00 0.02 0.05 0.07 0.10 0.12 TIME (XlO') co. (*(D I,.00oo 0.02 0.05 0.07 O.1 0.12 TIME (X101) "'t - -~TIE X1'

10. PROGRAM LISTING 45

C MARCH 24, 1986 C ==-===FOURTH ORDER RUNGE-KUTTA SCHEME —====== C C THIS PROGRAM GENERATES A COMPUTATIONAL MODEL OF THE C AUTOMOTIVE CLUTCH SYSTEM. IT IS MEANT TO REDUCE PRESENT C "TRIAL AND ERROR" METHODS USED IN CLUTCH DESIGN WHICH C 1HAVE NOT ONLY PROVED TO BE EXPENSIVE AND TIME CONSUMING, C BUT DO NOT ALLOW THE ENGINEER TO REALIZE THE POTENTIAL FOR C THE IMPROVED AND OPTIMAL DESIGNS WHICH COULD BE ACHIEVED. C C THE STRUCTURE OF THE PROGRAM CODING IS QUITE MODULAR, C ENABLING EASE OF MODIFICATION, ADDITION AND DELETION OF C VARIOUS SYSTEM PARAMETERS. IT ALSO INCORPORATES A FRIENDLY C ENVIRONMENT WITH 7hTE USER, EQUIPPED WITH MENUS AND EXAMPLES, C THEREFORE REQUIRING LITTLE KNOWLEDGE OF THE PROGRAM C TO RUN IT. THE USER CAN bAIN VALUABLE UNDERSTANDING OF C SYSTEM VARIABLES BY GENERATING A GRAPH VERSUS TIME, OR A C TABLE LISTING. ALSO, SYSTEM PARAMETERS CAN INTERRACTIVELY BE C MODIFIED BY THE USER AND THEN RERUN THE INTEGRATION. C THE NEW RESULTS WILL ONCE AGAIN BE AVAILABLE FOR POSTC PROCESSING. C C SUBROUTINE DESCRIPTIONS C DEFALT SETS DEFAULT VARIABLE CONSTANTS AS WELL AS C DEFAULTED COEFFICIENTS FOR PARAMETER EQUATIONS. C GRAPH GENERATES A X-Y GRAPH OF ANY TIME VARIANT C PARAMETER VERSUS TIME, AS SPECIFIED BY THE C USER C COUNTC COUNTS THE NUMBER OF CHARACTERS IN A PARTICULAR C CHARACTER STRING. CALLED BY GPAPH. C FILE ASSIGNS A FORTRAN UNIT NUMBER 10 A FILE NAME C WHICH IS SPECIFIED BY USER. CALLED BY PRINTA. C MODIFY ALLOWS USER TO MODIFY SYSTEM CONSTANTS, OR C SYSTEM PARAMETER EQUATION COEFFICIENTS BY C ENTERING DISCRETE DATA POINTS. CALLED BY MAIN. C PRINTA ALLOWS USER TO HAVE A TABLE OF THE TIME C DEPENDENT VARIABLES EITHER PRINTED AT THE C TERMINAL OR ONTO A SPECIFIED FILE. C C INTEGER CHR.YYQQSSMM,RR,GGPPICODE(14) REAL Q(12),M2,M4,11,I2,I3,14,IE,N3 REAL K2,K3,KHKDKE C COMMON/STATE/CHR COMMON/VEL/VRNEWVROLD COMMON /PAR/ M2,M4,11,12,13,14,IE., $TtiKCTHK5,THK4,THIK3,D2,D3,D5,RHUS.N3,UFTR. SRF.RC,A2.A1,AA, ICDUNT.DT.ddd COMMON /VARS/ TEFINPUT,K2,K3,KE,C2,C3,CEFDRDIFL C DATA,JQ / 0 / DATA YYQQ,GGMM.PP.RR.SS /'Y,'?'11,lG',#M''P','R','S' I C 1002 FORMAT (1/,' NEXT STEP? (ENTER "7" FOR MENU)') 1004 FORMAT (/,' AVAILABLE OPTIONS ARE::',//.' G GRAPH IWAOlTADI IC %IC T1MUC I/' M MnnllFY rCnNcNTANTS' _/

&' P PRINT VARIABLE TABLE',/,' R RUN INTEGRATION'./, &' S STOP') 1006 FORMAT (//.' ENTER DESIRED NUMBER OF TIME STEPS AND INCREMENT') 1010 FORMAT (//,' *+ INVALID OPTION -- TRY AGAIN **') 1020 FORMAT (//,' ** OUTPUT.OF QS FOR EACH TIME STEP CAN BE', &' FOUND ON',/.30A1,'**',///) 1022 FORMAT (//,' *** POSTPROCESSING OPTIONS ARE' &'NOT AVAILABLE UNTIL',/,' INTEGRATION', &' TECHNIQUE IS RUN (ENTER OPTION R) ***') 2000 FORMAT (Al) 2002 FORMAT (!5,F7.4) 2004 FORMAT (1X.I4) 20 FORMAT(2F9.4) 60 FORMAT(15) C C C C IFL = 0 CALL DEFALT C C C —-INITIALIZE NUMBER OF INDEP VARIABLES AND TIME —C N - 10 C C —-INITIALIZE INDEPENDENT VARABLES —C 0 0(1) = -10.5 0(2) = 0.0 0(3) - 0.0 0(4) = 0.0 o(5) - 0.0 0(6) = 0.0 0(7) - 0.0 0(8) - 52.36 0(9) = 0.0 0(10) = 0.0 Q(11)=Q(7)-Q(9) Q(12)=Q(8)-Q(10) C —-INITIALIZE COUNTER ICOUNT = 0 C C C —-PROCESS NEXT OPTION —C 105 WRITE (06, 1002) 1iO READ (05.2000) IANS C C —-LIST OPTIONS —C IF (IANS.NE. QQ) GO TO 115 WRITE (06. 1004) GO TO 105 C C —-STOP —C

115 IF (IANS.NE. SS) GO TO 120 GO TO 500 C C —-GRAPH —C 120 IF (IANS.NE. GG) GO TO 130 IF (IFL.NE. O) GO TO 122 WRITE (06,1022) GO TO 105 122 CALL GRAPH GO TO 105 C C —-MODIFY VARIABLE VALUES —C 130 IF (IANS.NE. MM) GO TO 140 IF (IFL.NE. O) GO TO 133 WRITE (06, 1022) GO TO 105 133 CALL MODIFY C C —-PRINT TABLE —C 140, IF (IANS.NE. PP) GO TO 150 IF (IFL.NE. 0) GO TO 144 WRITE (06,1022) GO TO 105 144 CALL PRINTA (Q,N) GO TO 105 C C —-RUN INTEGRATION TECHNIQUE —C 150 IF (IANS.NE. RR) GO TO 155 IFL = 1 REWIND 8 C C —-WRITE NUMBER OF TIME STEPS TO PLOT FILE —C WRITE (06,1006) READ (05,2002) NSTEP,TINC WRITE (06,200()) NSTEP,TINC WRITE (08,2004) NSTEP TFINAL=TINC*NSTEP WRITE (07,2019) TFINAL 2019 FORMAT (E10.4) C C —-FOR EACH TIME STEP —C WRITE (06,23) 23 FORMAT ('',' PRINT OUTPUT AT EVERY HOW MANY STEPS?') READ (05,22) ISTEP 22 FORMAT(12) C C STORE THE NUMBER OF PRINTED TIME STEP IN UNIT 07. NPSTP=NSTEP/ISTEP WRITE (07,2004) NPSTP C NS=O TAU=O.O C 79 DO 77 J=1,ISTEP

NS=NS+1 CALL KUTTA (Q,QDOT,TINC) TAU=TAU+TINC 77 CONTINUE C —-WRITE TIME, AND Q'S TO PLOT FILE —C WRITE (08,53) CHR,TAU,(Q(I), I=1,10) 53 FORMAT (2X,'STATE=',I4,10X,'TIME=',E10.4,/, $ 2X'( 1)=',EO.4,5X,'0(2)=',E10.4,5X,'Q(3)=',EIO.4,/,2X,'Q(4)', $ EIO.4,5X,'Q(5)=',E10.4,SX,'(6)=',E10.4,/,2X,'(7)=',EIO.4,5X, $'Q(8)=',EO.45X,'Q(9)=',EIO.4,/,2X,'Q( 10)=',E1.4,/) C Q( 11 )=O(7)-Q(9) Q( 12) Q(8)-Q( 10) C STORE THESE VALUES IN FILE PLOT.DATA FOR PLOT USAGE. WRITE (07,71) TAU,(Q(I),I=1,12) 71 FORMAT (13E15.7) C C IF (NS.GE.NSTEP) GO TO 105 GO TO 79 C C —-INVALID OPTION —C 155 WRITE (06,1010) GO TO 110 C C 500 WRITE (06, 1020) C STOP END C C SUBROUTINE DEFALT C C REAL M2,M4, I1,I2,I3,I4,IE,N3 REAL K2,K3,KH,KD,KE C COMMON /PAR/ M2,M4, 1,I2,I3,I4, IE, $THKC,THK5S,THK4,THK3,D2,D3,D5,RH,US,N3,UF,TR, $RE,RC,A2,A1,AA,ICOUNT,DT,dJJ COMMON /VARS/ TE,FINPUT,K2,K3,KE,C2,C3,CE,FD,RD, IFL C C —-SET VARIABLE CONSTANTS —C FD=O. O RD300. O A2=.281 A=1 i.O M2=10.0 M41.4 I1=169.48 12=0.5 I3=5.99 14=30.0 IE=200.0

THKS=3.8 THK4 =0. 6 THK3=3.8 D5= 1.0 RH=25.8 RS 13.49 UH=O. O N3= 100. O UF =0. 3 TR= 1. 000 RE=304.8 DT= 1.8 RC=93.75 C C RETURN END C C SUBROUTINE KUTTA (QQDOT,TINC) C INTEGER N,CHR REAL Q(20),QDOT(20),M2,M4,I1,I2,I3,.4.IEN3,0MEGA REAL K2,K3,KH,KD,KE,NN1 REAL Q1(10).Q2(10),Q3(10),DQI(10)'.DQ2(10),DQ3(10),DQ4(10).DQ(10) C COMMON /VARS/ TE,FINPUT,K2,K3,KEC2,C3,CE,FD,RD.IFL COMMON/STATE/CHR COMMON/VEL/VRNEW, VROLD COMMON /PAR/ M2,M4,I1,12,I3,I4,IE, $THKC.THK5S,THK4THK3.D2,D3.D5.RH.US,N3,UF,TR. $RE,RC,A2,A1,AA,ICOUNT,DT.JJU C C OMEGA=Q( 8 )-Q(6) C CALL TOROIN (UIMEGA) CALL FORCIN (T) C C —-STATE 1 TEST —C IF(Q(1).GF ( rT/A2)) GOTO 10 CALL STATE1 (U,(jI')OT) CHR= 1 C DO 700 I=1,10 DQ1( I )=TINC*QDOT( I ) Q1( )=Q( I )+DQi(I )/2.0 700 CONT I NUE C CALL STATE1 (QI,QDOT) C DO 710 I=1,10 DQ2( )=TINC*QDOT( I ) o2( I )oQ(I )+DQ2( I )/2.0 7 10 CONT I NUE C CALL STATE1 (0Q 2ODOT)

C DO 720 1=1,10 D3( I )=TINC*QDOT(I ) 03( I )(I)+D3( I ) 720 CONT I NUE C CALL STATE1 (Q3,QDOT) C DO 730 I=1,10 D04( I )-TINC*DOT( I ) DQ( I )=(DO1 (I )+2.0*(oQ2( I ) +DQ3(I ) )+DQ4( I ) )/6. Q(I )=Q(I)+DQ(I) 730 CONT I NUE RETURN C C —-STATE 2 TEST —C 10 IF(Q(3).GE..999) GO TO 20 CALL STATE2 (Q,QDOT) CHR=2 C DO 800 I=1,10 DQ1( I )=TINC*ODOT(I ) 01(1 )=Q( I )+DQI(I )/2.0 800 CONTINUE C CALL STATE2 (01,QDOT) C DO 805 1=1,10 DQ2(1 )=TINC*QDOT( I ) 02( I )Q(I )+DQ2( I )/2.0 805 CONT I NUE C CALL STATE2 (Q2,QD0T) C DO 810 -1, 10 DQ3( I )=TINC*QDOT(I ) 03(I )=(I)+D003(I ) 810 CONTINUE C CALL STATE2 (03,QDOT) C DO 815 1=1,10 DQ4(1 )=TINC*DOT( I ) DQ( I )(DQ1 (I )+2.0*(Do2( I )+D3( I ) )+DQ4( I ) )/6.0 Q( )=Q(I)+Do(I) 815 CONTINUE RETURN C C —-STATE 3 TEST —C 20 (4 )=0.0 IF (O(1).GE.-.001) GO TO 60 I COUNT= I COUNT+ i IF (ICOUNT.EQ. ) VRNEW=Q(8)-Q(10) IF (ICOUNT.EQ.1) ddd=O IF (ICOUNT.EQ.1) GO TO 500 C C NOATION: (1) dJJ=0, SLIPPING. C (3) JJJ=2, THE IMMEDIATE HISTORY IS ONE

C OF NON-SLIP. THE INTERFACE SURFACES C ARE BEGINNING TO SLIP. C (2) dJJ=1,. STICKING. C IF (JdJ.EQ.1) GO TO 200 IF ((VRNEW*VROLD).LT. 0.0) GO TO 200 C C 500 CALL STATE3 (QQDOT) IF (JJJ.EQ.O) CHR=3 IF (.JJd.EO.2) CHR-35 C DO 820 1=1,10 DQ1 ( I )=TINC*QDOT( I) Qo( I)=o( I )+DO1 (I)/2.0 820 CONT I NUE C CALL STATE3 (Q1,QDOT) C DO 825 1=1,10 DQ002(I )=TINCODOT( I ) 02( I )=Q( I )+D2(1 )/2.0 825 CONTINUE C CALL STATE3 (02,QDOT) C DO 830 1=1,10 DQ3(I )=TINC*ODOT( I ) 3( I )Q( I )+Q3( I) 830 CONT I NUE C CALL STATE3 (Q3,QDbT) C DO 835 I=1, 10 DQ4( I )=TINC*QDOT( I ) DQ( )=(DQ ( I )+2.0 (DQ2( I)+DQ3( I) )+DQ4( I) )/6.O Q(I)=Q(I)+DQ( I ) 835 CONTINUE C VROLD=VRNEW VRNEW=Q(8)-Q(10) C RETURN C C 200 TT1 = TE NNI = K3*((DT+A2*Q(1))-Q(3)) AA = I1*(TT1)/(I1+I2+I3+I4)-TT1 FT = NNI*UF*RC*2.0 IF(ABS(AA).GT. FT) dJJ=2 IF(ABS(AA).GT. FT) GO TO 500 IF(ABS(AA).LT. FT) JJJ=1 IF(dJJ. NE. 1) RETURN CALL STAT3C(Q.QDOT) VRNEW a 0.0 CHR = 33 c

DO 840 I=1,10O DOI(I)=TINC*ODOT(I) 1(I )=Q(I)+DQ( I )/2.0 840 CONTINUE C CALL STAT3C (Q1,QDOT) C DO 845 I-1,10 DQ2(X )TINC*QDOT(I) o2(1 )=Q( I )+Do2( 1)/2.0 845 CONTINUE C CALL STAT3C (Q2,QDOT) C DO 850 Izl,10 DQ3(I )=TINC*QDOT( I ) 03(I)=Q(I)+DQ3( I ) 850 CONTINUE C CALL STAT3C (Q3,QDOT) C DO 855 I=1.10 DQ4(I)=TINC*QDOT(I) DQ(I)=(DQI(I)+2.0*(DQ2(I)+DQ3(I))+DQ4(I))/6.0 ( I )=Q(I)+DQ(I) 855 CONTINUE RETURN C C C —-— BEGIN STATE4 FRICTION LOGIC —-- 60 0(2)=0.0 WRITE(8,9999) 0(8),Q(10) 9999 FORMAT(2X,'Q(8)=',EIO.4,5X,'Q(10)=',E10.4) WRITE(8,9998) VRNEW,VROLD 9998 FORMAT(2X,'VRNEW-',EO.4,4X,'VROLD=',E EIO.4) IF (JJJ.EQ.1) GO TO 205 IF((VRNEW*VROLD).LT. 0.0) GO TO 205 C C dJJ=O 505 CALL STATE4 (O,QDOT) IF (JJ. EQ.O) CHR =4 IF (JJd. EQ.2) CHR =45 C DO 860 1=1,10 DQ1(I)=TINC*QDOT(I) 1( I )=Q(I )+DQI(I)/2.0 860 CONTINUE C CALL STATE4 (Q1,QDOT) C DO 865 Il, 10 DQ2(1)=TINC*QDOT( I) Q2()=Q( I )+DQ2()/2.0 865 CONTINUE C CALL STATE4 (Q2,QDOT) C DO 870 I1=,10 DQ3(I)=TINC*QDOT( I.)

03(1 )=o(I)4DQ3U ) 870 CONTINUE C CALL STATE4 (Q3,QDOT) C DO 875 1=1,10 004()=TINCC*QDOT(I) DQ(I)=(DQ1(I)+2.O*(Q2(I )+DQ3(I) )+DQ4())/6.0 Q(I)=Q(I)+DQ(I) 875 CONTINUE C VROLO=VRNEW VRNEW=Q(8)-Q(10) C RETURN C C 205 TTI = TE NN1I K3*((DT+A2*Q(1))M0(3)) AA = I1*(TT1)/(I1+I2+I3+I4)-TT1 C C FT = NN1*UF*RC*2.0 C IF(ABS(AA).GT. FT) JJd=2 IF(ABS(AA).GT. FT) GO TO 505 IF(ABS(AA).LT. Fr) JJJ=1 IF(JJJ. NE. 1) RETURN C C CALL STAT4C (Q.QDOT) JJJ=i VRNEW c 0.0 CHR=43 C DO 880 1-1,10 DQ1(1)=TINCoQDOT(I) 01(1)=o( I)+oo( )/2.0 880 CONTINUE C CALL STAT4C (O1,ODOT) C 00 885 I=1,10 DQ2(1)=TINC*QOOT(I) Q2(I)=o(I)+D2( 1)/2.0 885 CONTINUE C CALL STAT4C (02.QDOT) C DO 890 1=1,10 003(1)=TINC*QDOT(I) 03(1)=Q(I))+DQ3(I 890 CON'TINUE C CALL STAT4C (Q3,.QOT) C DO 895 1-1,10 D04(1)=TINC*QDOT(I) DQ(I)=(DoQ( )+2.o* (oo2( DQ1)+DQ3(I))44())/6.0 o(I)=o(I)+oD(I)

895 CONTINUE RETURN END C C C = -== S 5= ===~ SUBROUTINE STATE1 (Q,QDOT) C Pressure plate not in contact with clutch: C clutch pedal starting return. REAL Q(10),QDOT(1O),K2,K3,KH,KD,KE,M2.M4,I1,I2,I3,I4,IE,N3 C COMMON /VARS/ TE,FINPUT,K2,K3,KE,C2,C3,CE,FD,RD,IFL COMMON /PAR/ M2,M4,1i,I2,I3,I4,IE, $THKC,THK5,THK4,THK3,D2,D3,D5,RH,US,N3,UF,TR, $RE,RC,A2,A1,AA,ICOUNT'DTJJJ C C BELSPR and ENGINE supply non-linear parameters K2 and C KE,CE respectively. C CALL BELSPR CALL ENGINE C C STATE #1 SET I - ROTATIONAL EQUATIONS OF MOTION C Flywheel/clutch cover/pressure plate/engine crank assembly C QDOT(7)=Q(8) QDOT(8)=TE/I1 C C STATE #1 SET II - ROTATIONAL EQUATIONS OF MOTION C Engine block C QODOT(5)=Q(6) QDOT(6)=-(RE/IE )*(CE*Q(6)+KE*Q() )+TE/IE C C STATE #1 SET I - TRANSLATIONAL EQUATIONS OF MOTION C Pressure Plate C THKT=9.34 QDOT(1)=Q(2) ODOT(2)=K2/(M2*A2)*(THKT- A2*Q( ))$FINPUT/(M2*A2)-(C2/M2)*Q(2) C QDOT(3) = 0. QDOT(4) = 0. QDOT(9) = 0. QDOT(10) = 0. RETURN END C C SUBROUTINE STATE2 (Q,QDOT) C Pressure plate In contact with clutch; flywheel C not in contact with clutch; slipping at clutch/pressure plate C interface assumed REAL Q(10),QDOT(10).K2,K3,KH,KD,KE,M2,M4,I1,I2.

& I3,I4.IE,N3,N1iK4 C COMMON /VARS/ TE,FINPUT,K2,K3,KE,C2.C3,CEFD, RD, II FL COMMON /PAR/ M2,M4,I1,2,I3,.T4,IE, $TtIKC.THKS, TtK4,THK3,D2,D3,D, RHUS,N3.UFTR, $RE,RC,A2,A1,AA.ICOUNT,DT,IJ C C BELSPR, MARCEL, and ENGINE supply non-linear C parameters (K2,C2), (K3,C3), and (KE,CE) C respectively. C CALL BELSPR CALL MARCEL CALL ENGINE C C N1, and T2 represent clutch normal C force, and clutch output torque respectively. C These quantities are used frequently in C subsequent equations and logic throughout the program. C Ni=K3*((DT+A2*Q(i))-Q(3))+C3*((A2*Q(2))-Q(4)) T1=NI*UF*RC C C STATE #2 SET I - ROTATIONAL EQUATIONS OF MOTION C Flywheel assembly C QDOT ( 7 )=( ) QDOT(8)=(TE-T1 )/I 1 C C STATE #2 SET II ROTATIONAL EQUATIONS OF MOTION C Outer clutch rotation C QDOT(9)=( 10) QDOT(l )=(T1)/(I2+I3+14) C C STATE #2 SET V - ROTATIONAL EQUATIONS OF MOTION C Engine block rotation. C QODOT(5)=Q(6) QDOT(s)=-(RE/IE )*(CE*Q(6)+KE*(5))+TE/IE C STATE #2 SET I - TRANSLATIONAL EQUATIONS OF MOTION C Pressure plate travel. C THKT=9.34 C QDOT(1)=Q(2) ODOT(2)=K2/(M2*A2)*(THKT-A2*Q(i))$FINPUT/(M2*A2)-(C2+C3)/M2*Q(2)+C3/(M2*A2)*Q(4)$K3/(M2*A2)*(DT+A2*Q(1)-Q(3)) C C STATE #2 SET II - TRANSLATIONAL EQUATIONS OF MOTION C Clutch displacement C ODOT(3)=Q(4) QDOT(4)=C3/M4*(Q(2)*A2-Q(4))+K3/M4*(DT+A2*Q(i)-Q(3)) RETURN END C C =3s= = s=-===========s=s=- s= =s= = z=-= = = = -= =-sX=- *=~==== -

SUBROUTINE STATE3 (Q,QDOT) C Jd-=O C Pressure plate and flywheel both in contact with clutch. C Interface surfaces are slipping. Clutch hub springs C are not fully compressed. C C JdJ=2 C The immediate history is one of non-slip. Therefore, C the clutch is contacting the flywheel and the pressure C plate. The Interface surfaces are BEGINNING to slip. C The direction of slip Is determined by the sign of C the quantity "AA" below. REAL Q(10),QDOT(iO),K2,K3,K4,KH,KD,KE,Ni $,M2,M4,I1, I2,I3,I4,IE,N3,N1MN2 C COMMON /VARS/ TE,FINPUT,K2,K3,KE,C2,C3,CE,FD,RD,IFL COMMON /PAR/ M2,M4,I1,12,I3,14,IE, $THKC,THK5,THK4,THK3,D2,D3,D5,RH,US,N3,UFTR, $RE,RC,A2,A1,AA,ICOUNT,DT,ddJ C BELSPR, MARCEL, and ENGINE supply non-linear C parameters (K2,C2), (K3,C3), and (KE,CE) C respectively. C CALL BELSPR CALL MARCEL CALL ENGINE C C N1, and T2 represent clutch normal C force, and clutch output torque respectively. C These quantities are used frequently in C subsequent equations and logic throughout the program. C N1i(K3*(DT+A2*Q(1)-Q(3))+C3*((A2*Q(2))-Q(4))) T2=NI*UF*RC*2.0 C C STATE #3 SET I - ROTATIONAL EQUATIONS OF MOTION C Equations check for direction of relative rotation (i.e. friction C torque direction) between flywheel assembly and clutch during C evaluation of flywheel assembly output values. C QDOT(7)=Q(8) IF (dJd.EQ.O) QDOT(8) = (TE-SIGN(1.O,Q(8)-Q(10))*T2)/I1 IF (JJJ.EQ.2) QDOT(8) = (TE+SIGN(1.O,AA)*T2)/I1 C C STATE #3 SET II - ROTATIONAL EQUATIONS OF MOTION C Equations check for direction of relative rotation C between flywheel assembly and clutch during C evaluation of clutch output characteristics. C QDOT(9)=Q(10) IF (JJd.EQ.0) QDOT(10)=(SIGN(1.O,Q(8)-Q(10))*T2)/(I2+13+I4) IF (JdJ.EQ.2) QDOT(10)=-(SIGN(i.O,AA)*T2)/(12+I3+14) C C STATE #3 SET V - ROTATIONAL EQUATIONS OF MOTION C Engine block rotation. C QDOT(5)=Q(6) QDOT(6)=-(RE/IE)*(CE*Q(6)+KE*Q(5))+TE/IE

C C STATE #3 SET I - TRANSLATIONAL EQUATIONS OF MOTION C Pressure plate C TIIKT9.34 C DOT ( ) Q(2) QDOT(2)=K2/(M2A2 ) * ( TIKT-A2*( i ) )$FINPUT/(M2*A2)-C2/M2Q(2(2)-C3/(M2*A2)*(Q(2)*A2-Q(4))SK3/(M2*A2)*(DT+A2*Q( )-Q(3)) C C Constraints for variables with fixed values; C specifically, the clutch is in contact with C the flywheel. C QDOT(3) = 0(4) QDOT(4)=O.O RETURN END C SUBROUTINE STAT3C (Q,QDOT) C Pressure plate and flywheel both in contact with clutch. C Interface surfaces are NOT slipping. Clutch hub springs C are NOT fully compressed; I.e. the stop pins are NOT C contacted. c==..=....= = =-=.=.=.. == = =.=.=.=.==.=======..... REAL o( iO),ODOT(O),K2,K3,K4,KH,KD,KE,N3,Ii,I2.I3,14,IE, & M2,M4,NI,N2 C COMMON /VARS/ TE,FINPUT,K2,K3,KE,C2,C3,CE,FD,RD,IFL COMMON /PAR/ M2,M4,11,I2,I3,I4,IE, $THKC,THKS,THK4,THK3.D2,D3,D5,RH.US,N3,UF.TR, $RE,RC,A2,AI,AA, COUNT,DT,JJ C C BELSPR, MARCEL, and ENGINE supply non-linear C parameters (K2,C2), (K3,C3), and (KE,CE) C respectively. C CALL BELSPR CALL MARCEL CALL ENGINE C C —-STATE #3C SET I - ROTATIONAL EQUATIONS OF MOTION C No slip is occurring at the clutch surface, therefore, the C flywheel assembly and clutch are treated as one rigid body. C For this reason, the respective rotational values are equated below. C QDOT(7)=Q(8) QDOT(8)=(TE)/(I 1+I2+I3+I4) oDOT( i0)=(TE)/(I1+12+I3+I4) ODOT(9)=Q( 10) C C STATE #3C SET IV - ROTATIONAL EQUATIONS OF MOTION C Engine block rotation. C QDOT(5)=Q(6) QDOT(6)=-(RE/IE)*(CE*Q(6)+KE*Q(5))+TE/IE

C STATE #3C SET I - TRANSLATIONAL EQUATIONS OF MOTION C Pressure plate C THKT=9.34 C ODOT(I)=0(2) OOoT(2)=(K2/(M2*A2))*(THKT-A2*Q(1))$FINPUT/(M2*A2)-C2/M2*Q(2)-.C3/(M2*A2)1Q(2)*A2-Q(4))$K3/(M2*A2)*(DT+A2*(i )-Q(3)) C C —-STATE #3C SET II - TRANSLATIONAL EQUATIONS —C C Constraints for variables with fixed values; C specifically, the clutch is in contact with C the flywheel. C QDOT(3)=Q(4) QDOT(4)=O.O RETURN END C SUBROUTINE STATE4 (QQOOT) C dJJ=.0 slipping. C dJJJ2. the immediate history is.-one of non-slip. C The interface surfaces are BEGINNING to slip. REAL Q(10),QDOT(10),K2,K3,K4,KHKD,KE,M2,M4,11,12, $ 13,14,IENi.N3 C COMMON /VARS/ TE,FINPUTK2,K3,KE,C2,C3,CEFDRDIFL COMMON /PAR/ M2,M4,11,12,13,14,IE, $THKCTHK5,THK4,THK3,D2,D3,05,RHUSN3,UF,TR, $RE,RC,A2,A1,AA,ICOUNT.DT,.JJ C CALL ENGINE CALL BELSPR CALL MARCEL Nl=(K3*(OT+A2*Q(1)-Q(3))+C3*((A2*Q(2))-Q(4))) T2=N1 *UF*RC*2.0 C C —-STATE #4 SET I - ROTATIONAL EQUATIONS —C QDOT(7)=Q(8) IF (JJJ.EQ.0) oDOT(8)=(TE-SIGN(1.o.o(8)-Q(10))*T2)/I1 IF (3JJ.EQ.2) QDOT(8)=(TE+SIGN(i.0,AA)*T2)/Ii C C C —-STATE #4 SET II & III - ROTATIONAL EQUATIONS —C QOoT(9)Q(i10) IF (JJJ.EQ.0) QDOT(iO)u(SIGN(i.00Q(8)-Q(10))*T2)/(12+13+14) IF (JJJ.EQ.2) QDOT(iO)m-(SIGN(1.0,AA)*T2)/(12+13+14) C C —-STATE #4 SET V - ROTATIONAL EQUATIONS —C QOOT(5)'Q(6) ODOT(6)=-(RE/IE)*(CEQ(6)+KE*Q(5) )-TE/IE C

C —-STATE #4 SET VI - TRANSLATIONAL EQUATIONS - CONSTRAINTS —C QDOT(1) = Q(2) QDOT(2) a 0.0 QDOT(3) = 0(4) QDOT(4) = 0.0 C RETURN END C SUBROUTINE STAT4C (Q,QDOT) C REAL Q(10).QDOT(10),KH,KD,KE.K2,K3, $ I1.I2,I3,I4,IEN3,M2,M4 C COMMON /VARS/ TE.FINPUT.K2,K3,KE,C2,C3,CE,FD.RD,IFL COMMON /PAR/ M2,M4,I,I2,13,14,IE, $THKC,THKS,THK4,THK3,D2,D3,DS,RH,US,N3,UF,TR, $RE,RC,A2,A1,AA,ICOUNT,DT,dJJ C CALL ENGINE C C QDOT(7)=Q(8) QDOT(8)(TE)/(I 1+I2+I3+I4) C CONSTRAINTS QDOT(9) =Q(10) QDOT( 10)=QOOT(8) QDOT( I )=Q(2) QDOT(2)-0.0 QDOT(3)Q(4 ) QDOT(4)=0.0 QDOT(5)=Q(6) ODOT(6) —(RE/IE)*(CE*Q(6)+KE*Q(5))+TE/IE RETURN END C SUBROUTINE TORQOIN (OMEGA) C C REAL TE,K2,K3.KHKD,KKE COMMON /VARS/ TE,FINPUT,K2,K3,KE,C2,C3,CE,FD,RD,IFL TE=10000.O RETURN END C SUBROUTINE FORCIN (T) C REAL T,FINPUT.K2,K3,KH,KD.KE COMMON /VARS/ TE.FINPUT.K2.K3.KE.C2.C3.CE.FD.RD.IFL

c FINPUT=O.O RETURN END C C SUBROUTINE BELSPR C REAL K2,C2,K3,KH,KD,KE COMMON /VARS/ TE,FINPUT,K2,K3,KE,C2,C3,CE,FD,RD, IFL C K2=200.0 C2=0.0 RETURN END C C SUBROUTINE MARCEL C = =-= = = = S-= -_ ===== == *== = =_ = C REAL K3,C3,K4,C4,K2,KH,KD,KE,M2,M4,I1,I2,I3,14,IE COMMON /VARS/ TE,FINPUT,K2,K3,KE,C2,C3,CE.FD,RD,IFL COMMON /PAR/ M2,M4,I1,I2,13,.4,IE, $THKC.THK5,THK4,THK3,D2,D3,D5.RHUS,N3,UF,TR, $RE,RC,A21,AA,ICOUNTDT,ddJ C K3=3000. 0 C3 = 0.0 RETURN END C C SUBROUTINE ENGINE C REAL KECE.K2,K3,KHKD COMMON /VARS/ TE,FINPUT,K2,K3,KE.C2,C3,CEFD,RD.IFL C KE=2000.0 CE=10.0 RETURN END C C SUBROUTINE GRAPH C LOGICAL*1 XT(50),YT(50),NAME(8) REAL X1(2500),YI(2500),DUM(12) INTEGER DD,TT,QQ,ICODE(12) C DATA DD,TT,QQ /'D','T','?' / DATA ICODE /'1','2','3','4','5','6','7', &'8','9'.'10''1 t''12'/ C

C 50 FORMAT (13E15.7) 70 FORMAT(5X,'ENTER X-AXIS LABEL') 80 FORMAT ( 50A 1 ) 90 FORMAT(SX,'ENTER Y-AXIS LABEL') C 1000 FORMAT (/,' WHICH VARIABLE? —ENTER "?" FOR MENU') 1002 FORMAT (/,' VARIABLES ARE::', & /.,' Qt',/,' 2 Q2',/,' 3 03', & /.' 4 O4',/,' 5 05',/.' 6 06', & /,' 7 07',/,' 7 Q7',/,' 8 08', & /,' 9 o9',/,' 10 otO',/,' 11 07-09' & /,' 12 Qs-Q10',/,' D DONE') 1006 FORMAT (//' *** INVALID ENTRY —TRY AGAIN ***') 2000 FORMAT (A2) 2002 FORMAT (1X,14) 2004 FORMAT (12) C C C C —-DEFINE VARIABLE TO BE PLOTTED VERSUS TIME —C 100 WRITE (06,1000) READ (05,2000) IANS C C —-PRINT LIST OF VARIABLES. IF REQUESTED —C IF (IANS.NE. QQ) GO TO 110 WRITE (06,1002) GO TO 100 C C —-IF USER DONE, RETURN —C 110 IF (IANS.EQ. DD) GO TO 500 C C —-TOROUE —C C IF (IANS.NE. TT) GO TO 10 C NN = 1 C GO TO 30 C C —-PROCESS ANY OF THE Q'S, IF CHOSEN —C 10 DO 20 I=1,12 IF (IANS.EQ. ICODE(I)) GO TO 25 20 CONTINUE WRITE (06,1006) GO TO 100 C 25 NN = I C C —-REWIND TEMPORARY FILE 7 —C 30 REWIND 7 C C —-READ FROM FILE 7 AND PACK X,Y ARRAYS —C READ (07,2019) TFINAL 2019 FORMAT (E10.4) RFAD (07.2002) N1

DO 40 I=1,N1 IF (NN.NE. 1) GO TO 35 READ (07,50) Xt(I),Yi(I),(DUM(d),J=l,11) YI(I)=-Y1(I) GO TO 40 C C —-READ TORQUE PLOT INFO —C FOR THIS VERSION, TORQUE IS NO LONGER A VARIABLE. C C READ (07,50) X1(I),Y1(I),(DUM(d),d=1,12),ICHR C GO TO 40 C 35 IF (NN.NE. 12) GO TO 37 READ (07,50) X1(I),(DUM(J),J=t, 1l),Y(I) GO TO 40 C 37 NZ = NN-1 NF = 12-NN READ (07,50) XI(I),(DUM(J),d=1,NZ),Y(I), (DUM(d).==I,NF) C 40 CONTINUE C C C —-GET AXIS LABEL INFORMATION —C WRITE(6,70) READ(5,80) XT CALL COUNTC(XT,50,NXT) WRITE(6,90) READ(5,80) YT CALL COUNTC(YT,40,NYT) C C XF=TFINAL/5. XM=O. O CALL PSCALE (3.,0.5,YM,YF,Y1,N1,1) CALL PLTOFS (XM.XF,YM,YF,1.,1.) CALL PLINE (X1.Y1,N,1,,,O0,1) CALL PAXIS (1.,1.,XT,-NXT,5.,0.,XM.XF,.5) CALL PAXIS (1.,1.,YT,NYT,3.,90.,YM.YF,.5) CALL PLTEND GO TO 100 500 CONTINUE C RETURN END C C SUBROUTINE COUNTC (STRN,MAXC,NCHR) C LOGICAL EQUC LOGICAL*1 STRN(1),BLNK DATA BLNK/''/ DO 20 I=t,MAXC IF(EQUC(STRN(I),BLNK)) GO TO 10 GO TO 20 10 IF(EQUC(STRN(I-1),BLNK)) GO TO 15

GO TO 20 15 NCHR=I-2 GO TO 30 20 CONTINUE NCHR I IF(EQUC(STRN( ),BLNK)) NCHR=I-1 30 RETURN END C C SUBROUTINE FILE (IUNIT,NAME) C C THIS SUBROUTINE ASSIGNS FORTRAN UNIT NUMBER, IUNIT, TO A FILE OF C WHICH IS CALLED NAME C LOGICAL*1 NAME(30) CALL SETLIO (IUNIT,NAME,&85) RETIURN 85 WRITE (06,87) 87 FORMAT (' ILLEGAL I/O UNIT OR FILE NAME') STOP END C C SUBROUTINE MODIFY C INTEGER ICODEA(13),ICODEB(f6),QQ REAL K2,K3,KH,KD,KE,M2,M4,I112, 13,1I4,IEN3 C COMMON /VARS/ TE,FINPUT,K2,K3,KE,C2,C3,CE,FD,RD, IFL COMMON /PAR/ M2,M4,I,12,13,14.IE, STHKC,THKS,THK4,THK3,D2,D3,D5,RH,US,N3,UFTR, $RE.RC,A2,A1,AA.ICOUNT,DT,JJJ C DATA dO / 0 / DATA ICODEA /'TE','F','K2','K3','KH','KD','KE'.'C2', &'C3','CE' E D'' B'.FD''RD' / DATA ICODEB'M2','M4','I1''I2','13','14','IE','RH', &'UH','US','N3','UF','TR','RE','RC','A2' / DATA DD /'0' / DATA OQ /'?' / C C 1000 FORMAT (/,' ON ONE INPUT LINE, ENTER VARIABLE CODE', &' AND MODIFY ID (C...CONSTANT; L...LEAST SQUARES FIT). & /,'TERMINATE INPUTS WITH "D <CR>". FOR', &' VARIABLE LISTING AND EXAMPLE ENTER "?"') 1002 FORMAT (//,' VARIABLES THAT MAY BE CHANGED ARE::', &' TE -- ENGINE TORQUE', & /,' FI -- CLUTCH PEDAL FORCE', & /,' K2 -- BELLEVILLE SPRING RATE', & /,' 3 -- CUSHION SPRING RATE', &.' KE -- ENGINE MOUNT STIFFNESS', & /,' H -- TORSIONAL CLUTCH SPRING RATE', & /,' -- TORSIONAL DRIVE TRAIN STIFFNESS', & /,' C2 -- BELLEVILLE LINKAGE DAMPING COEFFICIENT',

& /,' C3 -- CUSHION DAMPING COEFFICIENT', & /,' CE -- ENGINE MOUNT DAMPING COEFFICIENT', & /,' BD -- DRIVE TRAIN VISCOUS DAMPING COEFFICIENT'. & /,' FD -- ROLLING DRAG FORCE ON TIRE', & /,' RD -- TIRE ROLLING RADIUS', & //,' EXAMPLE INPUT SEQUENCE::' & /,' K3,C <CR> OR D <CR>'./) 1010 FORMAT (//,' ** INVALID VARIABLE CODE INPUT —', &'RE-ENTER LINE **') 2000 FORMAT (A2,A2) C C C 110 WRITE (06,1000) 120 READ (05,2000) IVAR,ICD IF (IVAR.NE. OQ) GO TO 200 WRITE (06, 1002) 00 = dQ+l IF (JO.GT. 1) GO TO 120 GO TO 110 C C —-CHECK IVAR FOR VALID INPUT —C 200 IC = 1 DO 300 Iw1,13 IF (IVAR.EQ. ICODEA(I)) GO TO 600 300 CONTINUE C 400 IC = 2 DO 500 I-1,16 IF (IVAR.EQ. ICODEB(I)) GO TO 600 500 CONTINUE WRITE (06,1010) GO TO 120 C C 600 IV = I CALL MODCON (IC,IV) RETURN END C C SUBROUTINE PRINTA (Q,N) C LOGICAL*1 NAME(30) REAL Q(N),M2,M4,N3,I 1,12,13.14, IE INTEGER TT,FF C COMMON /PAR/ M2,M4,I1,I2,13,14,IE, $THKC.THK5S,THK4,THK3,D2,D3,D5,RH,US,N3,UF,TR, $RE,RC,A2,A1,AA,ICOUNT,DT,JJJ C DATA TT,FF /'T','F' / 1000 FORMAT (/,' WOULD YOU LIKE YOUR TABLE TO BE PRINTED', &'-AT THIS TERMINAL',/,' OR SENT TO AN OUTPUT FILE', &' (T OR F)?') 1002 FORMAT (/,' PLEASE ENTER FILE NAME TO EMPTY DATA TABLE') 1004 FORMAT (//,' Q',G6X,'Q2',6X,'Q3',6X,'Q4',6X,'Q5',6X,'Q6',

& 6X,'Q7',6X,'Q8',6X,'Q9',5X,'QiO',/,' QI1',5X,'Q12', & 5X,'Q13',5X,'Q14',5X.'QDi',SX,'QD2',SX,'r,3', & 5X,'QD4'.5X,'QD5',5X.'QDG',/.' QD7',5X.'QD8',SX, &'Qo9',4X,'QDIO',4X,'QOD1',4X.'QDI2',4X,'QD13', & 4X.'QDi4',4X,'STATE') 2000 FORMAT (Ai) 2002 FORMAT (30A1) 2004 FORMAT (2X,2E10.4,/,2X,6El0.4/,/2X,6E10.4,/, $2X,2EtO.4,/,2X,I8) 2006 FORMAT (1X,I4) 2008 FORMAT (2X,7F10.4,/,2X,4F10.4,/2X,5F10.4,/,2X,18) C WRITE (06,1000) READ (05,2000) IANS C IF (IANS.EQ. TT) GO TO 200 IFILE = 4 WRITE (06,1002) READ (05,2002) NAME CALL FILE (IFILE,NAME) GO TO 300 C 200 IFILE - 6 C C C —-WRITE OUT TITLES FOR 0 AND QDOT OUTPUT —C —-ON DATA CHART FILE —C 300 WRITE (IFILE,1004) C C —-REWIND FILE 8 FOR READING —C REWIND 8 C C —-READ NUMBER OF STEPS —C READ (08,2006) NSTEP C C —-FOR EACH TIME STEP —C DO 500 I=1,NSTEP READ (08,2004) T,TOR,(Q(J),J=1, 10).IC1R WRITE (IFILE,2008) T,TOR,(Q(J).J=1, 10),ICHR 500 CONTINUE C RETURN END C C SUBROUTINE MODCON (IC,IV) C REAL RA(13)',RBI(7),RD(7),RB2(8),DUM(3) C COMMON /VARS/ RA,IDUM COMMON /PAR/ RB1,RD,RB2,DUM C 1000 FORMAT (/,' ENTER NEW CONSTANT VALUE') 1002 FORMAT (/.' *** INVALID IC IN MODCON *****')

2000 FORMAT (F9.4) C 300 IF (IC.NE. 1) GO TO 500 WRITE (06,1000) READ (05,2000) RA(IV) WRITE (06, 111) IV,RA(IV) ill FORMAT (' — VARS ARRAY ——',16,F10.4) GO TO 900 C 500 IF (IC.NE. 2) GO TO 800 WRITE (06, 1000) IF (IV.GT. 7) GO TO 700 READ (05,2000) RBI(IV) WRITE (06,222) IV,RBi(IV) 222 FORMAT (' — PAR ARRAY —-',I6,F10.4) GO TO 900 C 700 IN a IV-7 READ (05,2000) RB2(IN) WRITE (06,222) IN,RB2(IN) C 800 WRITE (06,1002) C 900 CONTINUE RETURN END

C MARCH 21, 1986 C = — c=FOURTH ORDER RUNGE-KUTTA SCHIEME':==='=-=== C C =~===J= STEP TORQUE C C THIS PROGRAM GENERATES A COMPUTATIONAL MODEL OF THE C AUTOMOTIVE CLUTCH SYSTEM. IT IS MEANI TO REDUCE PRESENT C "TRIAL AND ERROR" METHODS USED IN CLUTCH DESIGN WHICH C HAVE NOT ONLY PROVED TO BE EXPENSIVE AND TIME CONSUMING, C BUT DO NOT ALLOW THE ENGINEER TO REALIZE THE POTENTIAL FOR C THE IMPROVED AND OPTIMAL DESIGNS WHICH COULD BE ACHIEVED. C C THE STRUCTURE OF THE PROGRAM CODING IS QUITE MODULAR, C ENABLING EASE OF MODIFICATION, ADDITION AND DELETION OF C VARIOUS SYSTEM PARAMETERS. IT ALSO INCORPORATES A FRIENDLY C ENVIRONMENT WITH THE USER, EQUIPPED WITH MENUS AND EXAMPLES, C THEREFORE REQUIRING LITTLE KNOWLEDGE OF THE PROGRAM C TO RUN IT. THE USER CAN GAIN VALUABLE UNDERSTANDING OF C SYSTEM VARIABLES BY GENERATING A GRAPH VERSUS TIME, OR A C TABLE LISTING. ALSO, SYSTEM PARAMETERS CAN INTERRACTIVELY BE C MODIFIED BY THE USER AND THEN RERUN THE INTEGRATION. C THE NEW RESULTS WILL ONCE AGAIN BE AVAILABLE FOR POSTC PROCESSING. C C SUBROUTINE DESCRIPTIONS C DEFALT SETS DEFAULT VARIABLE CONSTANTS AS WELL AS C DEFAULTED COEFFICIENTS FOR PARAMETER EQUATIONS. C GRAPH GENFRATES A X-Y GRAPH OF ANY TIME VARIANT C PARAMETER VERSUS TIME, AS SPECIFIED BY THE C USER C COUNTC COUNTS THE NUMBER OF CHARACTERS IN A PARTICULAR C CHARACTER STRING. CALLED BY GRAPH. C FILE ASSIGNS A FORTRAN UNIT NUMBER TO A FILE NAME C WHICH IS SPECIFIED BY USER. CALLED BY PRINTA. C MODIFY ALLOWS USER TO MODIFY SYSTEM CONSTANTS, OR C SYSTEM PARAMETER EQUATION COEFFICIENTS BY C ENTERING DISCRETE DATA POINTS. CALLED BY MAIN. C PRINTA ALLOWS USER TO HAVE A TABLE OF THE TIME C DEPENDENT VARIABLES EITHER PRINTED AT THE C TERMINAL OR ONTO A SPECIFIED FILE. C C INTEGER CHR,YY.QQSS,MMRRGG.PP,ICODE(14) REAL Q(12),M2,M4,IlI2,I3,I4,IEN3 REAL K2.K3,KH,KDKE C COMMON/STATE/CHRNS,NTEATE COMMON/VEL/VRNEWVROLD COMMON /PAR/ M2,M4,11,12,13,14,IE, $THiKC.THK5,THK4,THK3,D2,D3,D5,RH,US,N3,UFTR, $RE,RC,A2.AiAAICOUNTDT,.JJ COMMON /VARS/ TEFINPUTK2,K3.KEC2,C3,CEFD,RDIFL C DATA JQ / 0 / DATA YY,QQ,GGMMPPRRSS /'Y','?','G','M',lP','R',lS' / C 1002 FORMAT (/1.' NEXT STEP? (ENTER "?" FOR MENU)')

1004 FORMAT (//,' AVAILABLE OPTIONS ARE::',//,' G GRAPH', &'VARIABLES VS TIME',/,' M MODIFY CONSTANTS',/, &' P PRINT VARIABLE TABLE'./.' R RUN INTEGRATION',/, &' S STOP') 1006 FORMAT (//,' ENTER DESIRED NUMBER OF TIME STEPS AND INCREMENT') 1010 FORMAT (//.' ** INVALID OPTION -- TRY AGAIN **') 1020 FORMAT (//,' ** OUTPUT OF QS FOR EACH TIME STEP CAN BE', &' FOUND ON',/,30A1,'**',///) 1022 FORMAT (//,' *** POSTPROCESSING OPTIONS ARE' &'NOT AVAILABLE UNTIL',/,' INTEGRATION', &' TECHNIQUE IS RUN (ENTER OPTION R) **') 2000 FORMAT (Al) 2002 FORMAT (15,F7.4) 2004 FORMAT (1X,I4) 20 FORMAT(2F9.4) 60 FORMAT(I5) C C C C IFL = 0 CALL DEFALT C C C —-INITIALIZE NUMBER OF INDEP VARIABLES AND TIME —C N 10O C C —-INITIALIZE INDEPENDENT VARABLES —C C C STEADY STATE VALUES FROM PREVIOUS RUN. C o(1) = 0.0 0(2) = 0.0 0(3) = 1.0 Q(4) = 0.0 0(5) = 0.01711 0(6) = -0.1987 Q(7) a 10.86 Q(8) - 52.76 Q(9) = 3.626 0(10) - 52.77 Q( 1 )=o(7)-0(9) ( 12)=Q(8)-o(10) C C C —-INITIALIZE COUNTER ICOUNT - 0 C C C —-PROCESS NEXT OPTION —C 105 WRITE (06,1002) 110 READ (05.,2000) IANS C C —-LIST OPTIONS —C IF (IANS.NE. QQ) GO TO 115

WRITE (06, 1004) GO TO 105 C C —-STOP —C 115 IF (IANS,NE. SS) GO TO 120 GO TO 500 C C- — GRAPH —C 120 IF (IANS.NE. GG) GO TO 130 IF (IFL.NE. O) GO TO 122 WRITE (06, 1022) GO TO 105 122 CALL GRAPH GO TO 105 C C —-MODIFY VARIABLE VALUES —C 130 IF (IANS.NE. MM) GO TO 140 IF (IFL.NE. O) GO TO 133 WRITE (06,1022) GO TO 105 133 CALL MODIFY C C —-PRINT TABLE —C 140 IF (IANS.NE. PP) GO TO 150 IF (IFL.NE. O) GO TO 144 WRITE (06,1022) GO TO tOS5 144 CALL PRINTA (Q,N) GO TO 105 C C —-RUN INTEGRATION TECHNIQUE —C 150 IF (IANS.NE. RR) GO TO 155 IFL = 1 REWIND 8 C C —-WRITE NUMBER OF TIME STEPS TO PLOT FILE —C WRITE (06,1006) READ (05,2002) NSTEP,TINC WRITE (06,2002) NSTEP,TINC WRITE (08,2004) NSTEP TFINAL=TINC*NSTEP WRITE (07,19) TFINAL 19 FORMAT(E10.4) C C —-FOR EACH tIME STEP —C WRITE (06,23) 23 FORMAT ('',' PRINT OUTPUT AT EVERY tlOW MANY STEPS 7') READ (05,22) ISTEP 22 FORMAT(I2) C C —-INPUT FORCE —----- WRITE (06,700) 700 FORMAT (''.' ENTER STEP TORQUE:')

READ (05,701) ATE 701 FORMAT (E10.4) WRITE (06,703) 703 FORMAT ('','STEP TORQUE AT WHAT TIME STEP?') READ (05,29) NTE 29 FORMAT (13) C C STORE THE NUMBER OF PRINTED TIME STEP IN UNIT 07. NPSTP=NSTEP/ISTEP WRITE (07.2004) NPSTP C NS=O TAU=O. O C 79 DO 77 d.1,ISTEP NS=NS+1 CALL KUTTA (Q,QDOT,TINC) U=TTAU+T INC 77 CONTINUE C —-WRITE TIME, AND Q'S TO PLOT FILE — C WRITE (08,53) CHR,TAU,(Q(I), I=1,10) 53 FORMAT (2X,'STATE=',4,10X,'TIME=',E10.4,/, $ 2X'Q( )=',E10.4,5X,'(2)-',EiO.4.5X.'(3)='. E1O.4,/,2X,'(4)=', $ EO.4.5X,'(5)-',EIO.4,5X,'(6)=',EO.4,/,2X,'(7)=',E1O.4,5X. $'(8)',EIO.4,5X,'(9)=',EIO.4,/,2X,'( 10)'.EO.4,/) C 0( 1 )=Q(7)-Q(9) Q( 12)=Q(8)-Q( 10) C STORE THESE VALUES IN FILE PLOT.DATA FOR PLOT USAGE. WRITE (07,71) TAU,(Q(I).I=1,12) 71 FORMAT (13E15.7) C C IF (NS.GE.NSTEP) GO TO 105 GO TO 79 C C —- INVALID OPTION —C 155 WRITE (06,1010) GO TO 110 C C 500 WRITE (06,1020) C STOP END C C SUBROUTINE DEFALT C C REAL M2,M4, 11,2,13,4, IE,N3 REAL K2,K3,KH,KD,KE C COMMON /PAR/ M2,M4.I1,I2,13,I4, IE, $THKC,THK5,THK4,THK3,D2,D3,DS,RH,US,N3, UF,TR, $RE,RC,A2,A1,AA,ICOUNT,DT.ddJ

COMMON /VARS/ TE,FINPUTK2,K3,KEC2,C3.CE,FDRDIFL C C —-SET VARIABLE CONSTANTS —C RD=300.0 A2a.281 Ai1i.0 M2=iO.O M22 10.0O M4= 1.4 112169.48 12=0.5 13=5.9S 14=30.0 IE=200.0 THK5=3.8 THK4O.6 THK3=3.8 D5=1.0 RH=25.8 RS13.49 UH=O. 0 N3= 100.0 UF=0.3 TRI.000 RE2304.8 DT=1.8 RC=93.75 C C RETURN END C C SUBROUTINE KUTTA (Q0.DOTTINC) C INTEGER N,CHR REAL Q(20),QDOT(20),M2,M4,11,12,13,14,IE.N3.OMEGA REAL K2,K3,KH.KDKE REAL QI(10),Q2(10).03(10),DQi(10),DQ2(10),DQ3(10),0Q4(10).DQ(i0) C COMMON /VARS/ TEFINPUTK2,K3.KE.C2,C3.CE.FD.RD.IFL COMMON/STATE/CHRR,NSNTE,ATE COMMON/VEL/VRNEW.VROLD COMMON /PAR/ M2,M41 1,12,13,14. IE, $THKC,THK5,THK4,THK3.D2,03,D5,RHUSN3,UFTR, $RE,RCA2,AI,AA,ICOUNTDT,JdJ C C OMEGAfQ(8)-0(6) C CALL TOROIN (OMEGA) IF (NS.GE.NTE) TE=TE+ATE C WRITE (06,11) TE C 11 FORMAT (.1','CURRENT TE IS:',EiO.4) C C —----— PROGRAM MODIFIED SO USER INPUT FORCE —------ crAl-l- FORCIN (T)

C C —-STATE I TEST —C IF(Q(1).GE.(-DT/A2)) GOTO 10 CALL STATEI (Q,QDOT) CHR= t C DO 700 I=1,10 DQ01( I )=TINC*QDOT( I ) Q1( I )Q( I )+DQI(I )/2.0 700 CONTINUE C CALL STATEI (Q1,QDOT) C DO 710 I=1,10 D002(I)=TINC*QDOT(I) 02(I)=Q(I )+DQ2(I )/2.0 710 CONTINUE C CALL STATE1 (Q2,QDOT) C DO 720 I=1,10 DQ3(I)=TINC*QDOT( I) Q3()=Q(I)+DQ3( I) 720 CONTINUE C CALL STATEI (03,QDOT) C DO 730 1=1,10 D04(1)=TINC*DOT( I) DQ(I )(D( I )+2.0*(DQ2(I )+DQ3( I) )+DQ4((I) ) )/6.0 O(I)=Q(I)+DQ(I) 730 CONTINUE RETURN C C —-STATE 2 TEST —C 10 IF(Q(3).GE..999) GO TO 20 CALL STATE2 (Q,QDOT) CHR=2 C DO 800 1=1,10 DQ1(I)=TINC*QDOT( I) Qi(l)=Q(I)+DQI(I)/2.0 800 CONTINUE C CALL STATE2 (Q1,QDOT) C DO 805 I=1,10 DQ2(I)=TINC*QDOT( I) 02(I )=Q(I )+DQ2(I )/2.0 805 CONTINUE C CALL STATE2 (Q2,QDOT) DO 810 I=l.10 DQ3(I=TINC*QDOT( I) 03(I)=Q(I)+DQ3(I) 810 CONTINUE C

CALL STATE2 (Q3.QDOT) C DO 815 I=1,10 DQ4( )rTINC*QOOT( ) DQ(I )-(DQo1( )+2.0*( DQ2(I)+DQ3( I) )+DQ4(I ))/G.0 O( I )Q( )+OQ( I ) 8 15 CONTINUE RETURN C C —-STATE 3 TEST —C 20 o(4 )=o. IF (Q(1).GE.-.001) GO TO 60 ICOUNT= I COUNT+ 1 IF (ICOUNT.EQ.i) VRNEW=Q(8)-Q(10) IF (ICOUNT.EQ.1) Jdd=0 C C NOATION: (1) dUd=0, SLIPPING. C (3) dJJ=2, THE IMMEDIATE HISTORY IS ONE C OF NON-SLIP. THE INTERFACE SURFACES C ARE BEGINNING TO SLIP. C (2) JJd=1, STICKING. C IF (UJd.EQ.1) GO TO 200 IF ((VRNEW*VROLD).LT. 0.0) GO TO 200 C C 500 CALL STATE3 (Q.QDOT) IF (dJU.EQ.O) CHR=3 IF (UJJ.EQ.2) CHR=35 C DO 820 I=1,10 DQi(I)=TINC*QDOT(I) Q( I )=( (I)+DO D(I )/2.0 820 CONTINUE C CALL STATE3 (01,QDOT) C DO 825 1=1,10 DQ2(I)=TINCkQDOT(I) Q2( I )=Q(I)+DQ2()/2.0 825 CONTINUE C CALL STATE3 (Q2,.DOT) C DO 830 1=1,10 DQ3( I)TINC*QDOT(I) Q3(I)=Q(I)+DQ3( I) 830 CONTINUE C CALL STATE3 (Q3,QDOT) C DO 835 I1,O10 D04( I )=TINC*QDOT( I ) DQ( I )=(DQo( I )2.0*(2()+.DQ3( I) )+DQ4( I) )/6.O O(I)Q(I )+DQ(I) 835 CONTINUE C VROLD=VRNEW VRNEW=Q(8)-Q(10)

C C C RETURN C C 200 TT1 = TE NNI = K3*((DT+A2*Q(1))-Q(3)) AA = I1*(TT1)/(I1+I2+I3+I4)-TT1 FT = NNI*UF*RC*2.0 C IF(ABS(AA).GT. FT) JJJ=2 IF(ABS(AA).GT. FT) GO TO 500 IF(ABS(AA).LT. FT) ddJ=1 IF(JdJ. NE. 1) RETURN C CALL STAT3C(Q,QDOT) VRNEW - 0.0 CHR w 33 C DO 840 I=1,10 DQ( I )=TINC*QDOT( I ) Q( I )=Q( I )+DQI(I )/2.0 840 CONTINUE C CALL STAT3C (Q1,QDOT) C DO 845 I=1,10 DO2( I )=TINC*QDOT( I ) 02(1)=Q(I)+DQ2( )/2.0 845 CONTINUE C CALL STAT3C (02,QDOT) C DO 850 I=1,10 DO3( I )=TINC*QDOT( I) Q3(I)sQ(I)+DQ3( I ) 850 CONTINUE C CALL STAT3C (Q3,QDOT) C DO 855 I=1,10 DQ4( I )-TINC*ODOT( I) DQ( I )=(DO1 ( I ).+2.0* (D2( I )+ (+D( I ) )Q4( I) )/6. Q(i )=Q(I )+DQ(I ) 855 CONTINUE RETURN C C C —-— BEGIN STATE4 FRICTION LOGIC —-- 60 O(2)=0.0 I COUNT = I COUNT+ 1 IF (ICOUNT.EQ. ) GO TO 205 WRITE(8,9999) O(8),Q(10) 9999 FORMAT(2X.'(8)',E10.4,5X,'(10)=',EIO0.4) IF (dJJ.EQ.2) WRITE(8,9998) VRNEW,VROLD 9998 FORMAT(2X,'VRNEW=', E 10.4,4X,'VROLD=', E 10O.4) IF (dJJ.EQ.1) GO TO 205 IF((VRNEW*VROLD).LT. 0.0) GO TO 205

C C JJJ=O 505 CALL STATE4 (Q,QDOT) IF (dJJJ.EQ.O) CIIR =4 IF (dJJ. EQ.2) CIIR,45 C DO 860 I-1, 10 D1 ( I )=TINC*QDOT( I) 01( I )=( I )+DQ1( I )/2.0 860 CONT I NUE C CALL STATE4 (Q1,QDOT) C a DO 865 It1, 10 DQ2( I )=TINC*QDOT( I ) o2( I )=Q(I )DQ2( I )/2.0 865 CONTINUE C CALL STATE4 (Q2,QDOT) C DO 870 I-, 10 DQ3( I )TINC*QDOT(I ) 03( )=Q( I )+D3( I ) 870 CONTINUE C CALL STATE4 (Q3,QDOT) C DO 875 1=1,10 DQ4( )=TINC*QDOT(I ) DQ( I )=(DQ( I )+2.0*(DO2( I )+DQ3( I) )+DQ4(I ) )/6.0 ( I )=Q( )+DQ( ) 875 CONTINUE C VROLD=VRNEW VRNEW-Q(8)-Q( 10) C C RETURN C 205 K3=3000.0 TTI = TE NN1 = K3*((DT+A2*Q(1))-Q(3)) AA = It1*(TT)/(I1+I2+I3+14)-TT1 C C FT = NNI*UF*RC*2.0 C IF(ABS(AA).GT. FT) Jdd=2 IF(ABS(AA).GT. FT) GO TO 505 IF(ABS(AA).LT. FT) dJJJu IF(JJJ. NE. 1) RETURN C C CALL STAT4C (Q,QDOT) JJJd= VRNEW = 0.0 CH1R=43 C DO 880 I1.10

DQI(I)=TINC*QDOT(I) 01(1)=Q(I )+Do0(I)/2. 880 CONTINUE C CALL STAT4C (QIQDOT) C 00 885 1=i,10 DQ2(1)=T1NC*QDoT(I) 02(1)=Q(I)+4DQ2()/2.0 885 CONTINUE C CALL STAT4C (02.QDOT) C 00 890 1=1,10 DQ3(I)=TINC*QDOT(I) 03(1)=Q(I)+0Q3(l) 890 CONTINUE C CALL STAT4C (03,QDOT) C DO 895 1=1,10 D04(1)=TINC*QDoT(I) DQ(I)=(OQi (I )+2.o*(D02(1I)+DQ3(I))+DQ4(I))/6.0 Q(I)=Q(I)+DQ(I) 895 CONTINUE RETURN END C C SUBROUTINE STATEI (Q,QOOT) C Pressure plate not in contact with clutch; C clutch pedal starting return. REAL (10),QODOT(10),K2K3,KHKDKE,M2,M4,.II,12,13,4,IEN3 C COMMON /VARS/ TE,FINPUT,K2,K3,KEC2,C3,CE,FDRDIFL COMMON /PAR/ M2,M4,11,12,.13;14,IE, $TNiKC,THK5,THK4,THK3,02,03,D5,RHUSN3,UF,TR, $RE,RCA2,A,AA,ICOUNT,DTJhJ C C BELSPR and ENGINE supply non-linear parameters K2 and C KECE respectively. C, CALL BELSPR CALL ENGINE C C STATE NI SET I - ROTATIONAL EQUATIONS OF MOTION C Flywheel/clutch cover/pressure plate/engine crank assembly C QDOT(7)=Q(8) QDOT(8)-TE/Ii C C STATE #1 SET II - ROTATIONAL EQUATIONS OF MOTION C Engine block C QDOT(5)=Q(6) QDOT(6)=-(RE/IE)*(CE*Q(6)+KEQ0(5))+TE/IE C

C Engine block rotation. C QDOT(5)=Q(6) QDOT(6)=-(RE/IE)*(CE*Q(6)KE*Q(5) )+TE/IE C C STATE #2 SET I - TRANSLATIONAL EQUATIONS OF MOTION C Pressure plate travel. C THKT=9.34 C QDOT(I)=Q(2) QDoT(2)=K2/(M2*A2)*(THKT-A2*Q(I))$FINPUT/(M2*A2)-(C2+C3)/M2*Q(2)+C3/(M2*A2)*Q(4)SK3/(M2*A2)*(DT4A2*Q(i)-Q(3)) C C STATE #2 SET II - TRANSLATIONAL EQUATIONS OF MOTION C Clutch displacement C QDOT(3)=Q(4) QDOT(4)=C3/M4*(Q(2)*A2-Q(4))+K3/M4*(DT+A2*Q(i)-Q(3)) RETURN END C SUBROUTINE STATE3 (QOQDOT) C JJJ~O C Pressure plate and flywheel both in contact with clutch. C Interface surfaces are slipping. Clutch hub springs C are not fully compressed. C C JJi32 C The Immediate history is one of non-slip. Therefore. C the clutch is contacting the flywheel and the pressure C plate. The interface surfaces are BEGINNING to slip. C The direction of slip is determined by the sign of C the quantity "AA" below. REAL Q(1O),QDOT(1O),K2.K3.K4,KH,KD,KENi $,M2.M4.I1,I2, I3,14, IEN3.N1MN2 C COMMON /VARS/ TE,FINPUTK2,K3,KEC2.C3,CEFDRD,IFL COMMON /PAR/ M2,M4.11,12,13,14.IE. $THKC.THKS,THK4,THK3,D2,D3,D5,RH.US,N3.UFTR, $RERC,A2.AI,AA,ICOUNT,DT.JJJ C BELSPR, MARCEL. and ENGINE supply non-linear C parameters (K2,C2), (K3,C3). and (KE.CE) C respectively. C CALL BELSPR CALL MARCEL' CALL ENGINE C C NI, and T2 represent clutch normal C force, and clutch output torque respectively. C These quantities are used frequently in C subsequent equations and logic throughout the program. C Ni=(K3*(DT+A2*Q(1)-Q(3))+C3*((A2*Q(2))-Q(4))) T2=N1*UF*RC*2.0

C STATE #1 SET I - TRANSLATIONAL EQUATIONS OF MOTION C Pressure Plate C THKT=9.34 QDOT( )=Q(2) QDOT(2)=K2/(M2A2)*(THKT-A2*Q( 1))$FINPUT/(M2*A2)-(C2/M2)*Q(2) C QDOT(3) = 0. QDOT(4) = O. QDOT(9) = O. QDOT(10) a O. RETURN END C C SUBROUTINE STATE2 (Q,QDOT) C Pressure plate in contact with clutch: flywheel C not in contact with clutch; slipping at clutch/pressure plate C interface assumed REAL Q(1O),QDOT(lO),K2,K3,KH,.KDKE,M2,M4, I1, I2, & I3.I4.IEN3,NI,K4 C COMMON /VARS/ TE,FINPUT,K2,K3,KE,C2,C3,CE,FD,RD. IFL COMMON /PAR/ M2,M4, I1,2,3,14, IE, STHKC,THK5,THK4,THK3,D2.13,D5,RH,US,N3,UF,TR, $RE,RC,A2,A1,AA,ICOUNT,DT,JJJ C C BELSPR, MARCEL, and ENGINE supply non-linear C parameters (K2,C2), (K3,C3), and (KE.CE) C respectively. C CALL BELSPR CALL MARCEL CALL ENGINE C C NI, and T2 represent clutch normal C force, and clutch output torque respectively. C These qnetrif, I,,' -, I- teFd frequlentlv In C subsequent equations and logic throughout the program. C NI=K3*((DT+A2*Q(i))-Q(3))+C3*((A2*Q(2))-Q(4)) TI=NI*UF*RC C C STATE #2 SET I - ROTATIONAL EQUATIONS OF MOTION C Flywheel assembly C QDOT(7)=Q(8) QDOT(8)-(TE-TI )/I C C STATE #2 SET 11I ROTATIONAL EQUATIONS OF MOTION C Outer clutch rotation C QDOT(9),Q( 10) QDOT(lO)=(TI)/(1I2+3+14) C C STATE #2 SET V - ROTATIONAL EQUATIONS OF MOTION

C C STATE #3 SET I - ROTATIONAL EQUATIONS OF MOTION C Equations check for direction of relative rotation (i.e. friction C torque direction) between flywheel assembly and clutch during C evaluation of flywheel assembly output values. C QDOT(7) Q(8) IF (JJJ.EQ.O) QDOT(8) = (TE-SIGN(1.0,Q(8)-Q(10))*T2)/I1 IF (dJd.EQ.2) QDOT(8) = (TE+SIGN(i.0.AA)*T2)/I1 C C STATE #3 SET II - ROTATIONAL EQUATIONS OF MOTION C Equations check for direction of relative rotation C between flywheel assembly and clutch during C evaluation of clutch output characteristics. C QDOT(9)=Q( 10) IF (dJJ.EQ.O) QDOT(1O) (SIGN(1.0,Q(8)-Q(10))*T2)/(I2+13+I4) IF (dJJ.EQ.2) QDOT(1O)=-(SIGN(.O,AA)*T2)/(I 2+I3+I4) C C STATE #3 SET V - ROTATIONAL EQUATIONS OF MOTION C Engine block rotation. C QDOT(5)=Q(6) QDOT(6)=-(RE/IE)*(CE*Q(6)+KE*Q(5))+TE/IE C C STATE #3 SET I - TRANSLATIONAL EQUATIONS OF MOTION C Pressure plate C THKT=9.34 C QDOT(1)=Q(2) QDOT(2)=K2/(M2*A2)*(THKT-A2*Q(1))$FINPUT/(M2*A2)-C2/M2*Q(2)-C3/(M2+A2)*(Q(2)*A2-Q(4))$K3/(M2*A2)*(DT+A2*Q(1)-0(3)) C C Constraints for variables with fixed values; C specifically, the clutch is in contact with C the flywheel. C QDOT(3) = 0(4) QDOT(4)=O.O RETURN END C SUBROUTINE STAT3C (O,QDOT) C Pressure plate and flywheel both in contact with clutch. C Interface surfaces are NOT slipping. Clutch hub springs C are NOT fully compressed; I.e. the stop pins are NOT C contacted. REAL Q(1O),QDOT(1O),K2,K3,K4.KH.KD.KE,N3,Ii.12. I3,14.IE, & M2,M4,N1,N2 C COMMON /VARS/ TE,FINPUT,K2,K3,KE,C2,C3,CE,FD,RD.IFL COMMON /PAR/ M2,M4,11,12,3,I14,IE, $THKC,THK5,THK4,TtHK3,D2,D3,D5,RH.US.N3,UF.TR, $RE,RC,A2,A1,AA,ICOUNTDT.JdJJ C

C BELSPR, MARCEL, and ENGINE supply non-linear C parameters (K2,C2), (K3,C3), and (KECE) C respectively. C CALL BELSPR CALL MARCEL CALL ENGINE C C —-STATE #3C SET I - ROTATIONAL EQUATIONS OF MOTION C No slip is occuring at the clutch surface, therefore, the C flywheel assembly and clutch are treated as one rigid body. C For this reason, the respective rotational values are equated below. C QDOT(7)-Q(8) QDOT(8)=(TE)/(1 +I2+I3+I4) QDOT(IO)=(7E)/(I lI2+[3+14) QDoT(9)=Q(10) C C STATE N3C SET IV - ROTATIONAL EQUATIONS OF MOTION C Engine block rotation. C QDOT(5)=Q(6) ODOT(6)=-(RE/IE)*(CE*Q(6)+KE*Q(5))-TE/IE C C STATE N3C SET I - TRANSLATIONAL EQUATIONS OF MOTION C Pressure plate C THKT=9.34 C QDOT( 1)=Q(2) QDOT(2)=(K2/(M2*A2) )*(THKT-A2*(1))$FINPUT/(M2*A2)-C2/M2*Q(2)-C2/(M2*A2)*(Q(2)*A2-Q(4))SK3/(M2*A2)*(DT+A2*Q( )-Q(3)) C C —-STATE #3C SET II - TRANSLATIONAL EQUATIONS —C C Constraints for variables with fixed values; C specifically, the clutch is in contact with C the flywheel. C QOOT(3)=Q(4) QUOT(4)=O.O RETURN END C SUBROUTINE STATE4 (Q,QDOT) C JJJ=0, slipping. C JJJ=2, the immediate history is one of non-slip. C The interface surfaces are BEGINNING to slip. REAL Q(10),QDOT(10),K2,K3,K4,KHKDKEM2,M4,I1,12. $ 13,14,IENIN3 C COMMON /VARS/ TEFINPUT,K2,K3,KE,C2,C3,CEFDRD,IFL COMMON /PAR/ M2,M4,11,12,13,14,YE, $THKC,THK5,THK4,THK3.02,03,D5,RH,USN3,UFTR, $RERCA2,A1.AAICOUNT.DTJJJ C

CALL ENGINE CALL BELSPR CALL MARCEL Nt=(K3*(DT+A2*Q( )-Q(3))+C3*((A2*(2))-Q(4)) ) T2=N1iUF*RC*2.0 C C —-STATE #4 SET I - ROTATIONAL EQUATIONS —C QDOT(7)=Q(8) IF (JJJ.EQ.O) QDOT(8)=(TE-SIGN(1.O,Q(8)-Q(10))*T2)/I1 IF (dJJ.EQ.2) QDOT(8)-(TE+SIGN(1.0,AA)*T2)/I1 C C C —-STATE #4 SET II & III - ROTATIONAL EQUATIONS —C QDOT(9)=Q(10) IF (dJJ.EQ.O) ODOT(10)=(SIGN(1.0,(8)-Q(1O))+T2)/(12+I3+I4) IF (JdJ.EQ.2) QDOT(1O)"-(SIGN(1.O,AA)*T2)/(12+I3+I4) C C —-STATE #4 SET V - ROTATIONAL EQUATIONS —C QDOT(5)-Q(6) QDOT(6) - (RE/IE)*(CE*Q(6)+KE*Q(5) )+TE/IE C C —-STATE #4 SET VI - TRANSLATIONAL EQUATIONS - CONSTRAINTS —C QDOT(1) - Q(2) QDOT(2) = 0.0 ODOT(3) a Q(4) QDOT(4) = 0.0 C RETURN END C SUBROUTINE STAT4C (Q,QDOT) C C REAL (10),QODOT(10),KH,KD.KE,K2,K3, $ I1,I2,I3,I4,IE,N3,M2,M4 C COMMON /VARS/ TE,FINPUT,K2,K3,KE,C2,C3,CE,FD,RD.IFL COMMON /PAR/ M2,M4,I1,12,I3,I4,IE, $THKC,THK5S,THK4, THK3,D2,D3,D5,RH,US,N3,UF.TR, $RE,RC,A2,A1,AA, ICOUNT,DT,JJJ C CALL ENGINE CALL MARCEL CALL BELSPR C C ODOT(7)=0(8) oDOT(8)=(TE)/(I 1+I2+I3+14) C CONSTRAINTS ODOT(9)o( 10) QDOT(10)=QDOT(8) ODOT( I )=Q(2) ODOT ( 2)=0.0

ODOT(3)=Q(4) QDOT(4)0.0 QDOT(5)=Q(6) QDOT(G)=-(RE/IE)*(CE*Q(6)+KE*Q(5))+TE/IE RETURN END SUBROUTINE TOROIN (OMEGA) C REAL TE,K2,K3,KH,KD,KE COMMON /VARS/ TE,FINPUT,K2,K3,KE,C2,C3,CE,FD,RD,IFL TE=10000.0 RETURN END C C SUBROUTINE FORCIN (T) C C=S======== ====S_=== = REAL T,FINPUT,K2,K3.KHKD,KE COMMON /VARS/ TE,FINPUT,K2,K3,KE,C2,C3,CE,FD,RD,IFL C FINPUT=O.O RETURN END SUBROUTINE BELSPR C REAL K2,C2,K3,KH,KD,KE COMMON /VARS/ TE,FINPUT,K2,K3,KE,C2,C3,CE,FD,RD,IFL C K2=200.0 C2=0.0 RETURN END C C SUBROUTINE MARCEL C 3=3=3sr= ===s-========_-=_=======33=3=3=s C REAL K3,C3,K4,C4,K2,KH,KD,KE,M2,M4,1l,I2,13,14,IE COMMON /VARS/ TE,FINPUT,K2,K3,KE,C2,C3,CE,FD,RD,IFL COMMON /PAR/ M2,M4,I1,12,13,14,IE, $THKC,THK5S,THK4,THK3,D2,D3,D5,RH,US,N3,UF,TR, $RE,RCA2,A1,AA,ICOUNTDTdJJ C K3=3000.0 C3 = 0.0 RETURN END C

C SUBROUTINE ENGINE C C REAL KECE,K2,K3.KtI,KD COMMON /VARS/ TE,FINPUTK2.K3,KEC2,C3,cE,FD,RD,IFL C KE x2000.0 CE=10.0 RETURN END C C SUBROUTINE GRAPH C LOGICALSi XT(50),YT(50),NAME(8) REAL XI(2500),YI(2500),DUM(12) INTEGER DD,TT,0QICODE(12) C DATA DDTTQQ /'D','T','?' / DATA IC3DE /'1','2','3 S,'4 lf5','6 i,'7 8,,,~~18 919 Jq v#1'q'lilq'12 C C 50 FORMAT (3F15 7) 70 FORMAT(5X,'ENTER X-AXIS LABEL') 80 FORMAT(SOA1) 90 FORMAT(5X,'ENTER Y-AXIS LABEL') C 1000 FORMAT (I,' WHICH VARIABLE? —ENTER "7" FOR MENU') 1002 FORMAT (U,' VARIABLES ARE::', & I / 02',/,' 3 03', & 1 4 4',/,' 5 05',!,' 6 06', & 7 07',/,' 7 Q7',/,' 8 08', & 9 09',!.' 10o 010',,' 11 07-09' & I, Q12 0-0',/,' D DONE') 1006 FORMAT (//' $ INVALID ENTRY —TRY AGAIN ***') 2000 FORMAT (A2) 2002 FORMAT (1X,14) 2004 FORMAT (12) C C C C —-OEFINE VARIABLE TO BE PLOTTED VERSUS TIME —C 100 WRITE (06,1000) READ (05,2000) IANS C C —-PRINT LIST OF VARIABLES, IF REQUESTED —C IF (IANS.NE. 00) GO TO 110 WRITE (06.1002) GO TO 100 C C —-IF USER DONE, RETURN —110 IF (IANS.EO. DD) GO TO 500

C C —-TORQUE —C C IF (IANS.NE. TT) GO TO 10 C NN a I C GO TO 30 C C —-PROCESS ANY OF THE Q'S. IF CHOSEN —C 10 00 20 I=1,12 IF (IANS.EQ. ICODE(I)) GO TO 25 20 CONTINUE WRITE (06.1006) GO TO 100 C 25 NN = I C C —-REWIND TEMPORARY FILE 7 —C 30 REWIND 7 C C —-READ FROM FILE 7 AND PACK X,Y ARRAYS —C READ (07,19) TFINAL 19 FORMAT (E10.4) READ (07.2002) Ni DO 40 Inl.NI IF (NN.NE. i) GO TO 35 READ (07,50) X1(I), Y(I).(DUM(J),d-i,1) YI(I)w-YI(I) GO TO 40 C C —- READ TORQUE PLOT INFO —C FOR THIS VERSION, TORQUE IS NO LONGER A VARIABLE. C C READ (07,50) Xi(I),.1(I),(DUM(J),=1, 12),ICHR C GO TO 40 C 35 IF (NN.NE. 12) GO TO 37 READ (07.50) X1(I),(DUM(J).J=1, 11),Y(I) GO TO 40 C 37 NZ = NN-1 NF a 12-NN READ (07,50) Xi(I),(DUM(dJ=1,Nz),v1(I),(DUM(UJ),MJ=NF) C 40 CONTINUE C C C —-GET AXIS LABEL INFORMATION —C WRITE(6.70) READ(5,80) XT CALL COUNTC(XT,50,NXT) WRITE(6,90) READ(5,80) VT CALL COUNTC(YT,40,NYT) C 22 FORMAT (FiO.4)

XF=TFINAL/5. XM=O.O C CALL PSCALE (3.,0.5,YMYF,VYNi,1) CALL PLTOFS (XMXF,YM,YF,1.,l.) CALL PLINE (X1,YINI,1.o.o,1) CALL PAXIS (1.,1.,XT,-NXT,,5..,.XM,XF,.5) CALL PAXIS (1.1.,YTNYT.3..90.YM.YF,.5) CALL PI-TEND GO TO 100 500 CONTINUE C RETURN END C C SUBROUTINE COUNTC (STRNMAXC,NCHR) C= — C LOGICAL EQUC LOGICAL*1 STRN(1),BLNK DATA BLNK/''/ DO 20 11I,MAXC IF(EQUC(STRN(I),BLNK)) GO TO 10 GO TO 20 10 IF(EQUC(STRN(I-1),BLNK)) GO TO 15 GO TO 20 15 NCHR=1-2 GO TO 30 20 CONTINUE NCHR=I IF(EQUC(STRN(I),BLNK)) NCHR=I-1 30 RETURN END C C SUBROUTINE FILE (JUNITNAME) C THIS SUBROUTINE ASSIGNS FORTRAN UNIT NUMBER, IUNIT, TO A FILE OF C WHICH IS CALLED NAME C LOGICAL4i NAME(30) CALL SETLIO (IUNITNAME,&85) RETURN 85 WRITE (06,87) 87 FORMAT (' ILLEGAL 1/0 UNIT DR FILE NAME') STOP END C C SUBROUTINE MODIFY C' == ======= uw== -— ======== = = === == = ===== ===================== INTEGER ICODEA(13),ICDDEB(i6),QQ REAL K2,K3,KHKD.KE,M2.M4, 1,12,I3,14,IE.N3 C

COMMON /VARS/ TE,FINPUT,K2,K3,KE,C2,C3,CE,FD,RD,IFL COMMON /PAR/ M2,M4, I1,I2,13,I4, IE, $THKC,THK5,THK4,THK3,D2,D3,D5,RH,US,N3,UF,TR, $RE,RC,A2,A,AA,ICOUNT,DT,JJJ C DATA JQ / 0 / DATA ICODEA /'TE','FI','K2''K3','KH','KD','KE','C2', &'C3','CE''BD''FD','RD' DATA ICODEB /'M2','M4''I'I'2','13','14','IE','RH', &'UH','US','N3','UF','TR','RE','RC','A2' / DATA DD /'D' / DATA QQ /'?' / C C 1000 FORMAT (/,' ON ONE INPUT LINE, ENTER VARIABLE CODE', &' AND MODIFY ID (C...CONSTANT; L...LEAST SQUARES FIT).', & /,'TERMINATE INPUTS WITH "D <CR>"., FOR', &' VARIABLE LISTING AND EXAMPLE ENTER "?"') 1002 FORMAT (//.' VARIABLES THAT MAY BE CHANGED ARE::', & /,' TE -- ENGINE TORQUE', & /,' FI -- CLUTCH PEDAL FORCE', & /,' K2 -- BELLEVILLE SPRING RATE', & /,' K3 -- CUSHION SPRING RATE', & /,' KE -- ENGINE MOUNT STIFFNESS', & /,' KH -- TORSIONAL CLUTCH SPRING RATE', & /,' KD -- TORSIONAL DRIVE TRAIN STIFFNESS', & /,' C2 -- BELLEVILLE LINKAGE DAMPING COEFFICIENT', & /,' C3 -- CUSHION DAMPING COEFFICIENT', & /,' CE -- ENGINE MOUNT DAMPING COEFFICIENT', & /,' BD -- DRIVE TRAIN VISCOUS DAMPING COEFFICIENT', & /,' FD -- ROLLING DRAG FORCE ON TIRE', & /,' RD -- TIRE ROLLING RADIUS', & //,' EXAMPLE INPUT SEQUENCE::' & /,' K3,C <CR> OR D <CR>',/) 1010 FORMAT (//,' ** INVALID VARIABLE CODE INPUT —', &'RE-ENTER LINE **') 2000 FORMAT (A2,A2) C C C 110 WRITE (06,1000) 120 READ (05,2000) IVAR,ICD IF (IVAR.NE. QQ) GO TO 200 WRITE (06, 1002) JOQ = d+1 IF (dQ.GT. 1) GO TO 120 GO TO 110 C C —-CHECK IVAR FOR VALID INPUT —C 200 IC = I DO 300 I-1,13 IF (IVAR.EQ. ICODEA(I)) GO TO 600 300 CONTINUE C 400 IC - 2 DO 500 It1,16 IF (IVAR.EQ. ICODEB(I)) GO TO 600 500 CONTINUE WRITE (06, 1010)

GO TO 120 C C 600 IV = I CALL MODCON (IC,IV) RETURN END C C SUBROUTINE PRINTA (Q.N) C LOGICAL*1 NAME(30) REAL Q(N),M2,M4,N3, I, I2,13,I4,IE INTEGER TT,FF C COMMON /PAR/ M2,M4,11,I2,I3,14,IE, $THKC,THK5,THK4,THK3,D2,D3.D5,RH,US,N3,UF,TR, $RE,RC,A2,AI,AA,ICOUNT,DT,dJJ C DATA TT,FF /'T','F' / 1000 FORMAT (/,' WOULD YOU LIKE YOUR TABLE TO BE PRINTED', &' AT THIS TERMINAL',/,' OR SENT TO AN OUTPUT FILE', &' (T OR F)?') 1002 FORMAT (/,' PLEASE ENTER FILE NAME TO EMPTY DATA TABLE') 1004 FORMAT (//,' Q1',6X,'Q2',6X,'Q3',6X,'Q4',6X,'Q5',6X,'Q6', & 6X,'Q7',6X,'8',X,'Q9',5X,'10',/,' Q11',5X,'Q12', & 5X,'Q13',5X,'Q14',5X.'QO1',5X,'QD2',5X,'QD3', & 5X,'QD4',5X, QD5'.5X'QD6',/.' QD7',5X,'QD8',5X, &'QD9',4X,'QDIO',4X,'QD11',4X.'QD12',4X,'QD13', & 4X,'QD14',4X,'STATE') 2000 FORMAT (Al) 2002 FORMAT (30A1) 2004 FORMAT (2X,2E10.4,/,2X,6E10.4,/,2X,6E10.4,/, $2X,2E10.4,/,2X,I8) 2006 FORMAT (1X,14) 2008 FORMAT (2X,7F10.4,/,2X,4F10.4,/2X.5F10.4,/,2X,18) C WRITE (06, 1000) READ (05,2000) IANS C IF (IANS.EQ. TT) GO TO 200 IFILE = 4 WRITE (06,1002) READ (05.2002) NAME CALL FILE (IFILE,NAME) GO TO 300 C 200 IFILE = 6 C C C —-WRITE OUT TITLES FOR 0 AND ODOT OUTPUT —C —-ON DATA CHART FILE —C 300 WRITE (IFILE,1004) C C —-REWIND FILE 8 FOR READING —C DFWIND a

C C —-READ NUMBER OF STEPS —C READ (08,2006) NSTEP C C —-FOR EACH TIME STEP —C DO 500 Ii1,NSTEP READ (08,2004) T,TOR,(Q(J),-1,1i0),ICHR WRITE (IFILE,2008) T,TOR,(Q(J),J1,I10),ICHR 500 CONTINUE C RETURN END C C SUBROUTINE MODCON (ICIV) C C S===u=====m3~~~=. D =.= = =.. = REAL RA( 13),RB1(7),RD(7),RB2(8),DUM(3) C COMMON /VARS/ RA,IDUM COMMON /PAR/ RBI,RD,RB2,DUM C 1000 FORMAT (/,' ENTER NEW CONSTANT VALUE') 1002 FORMAT (/,' *** INVALID IC IN MODCON *****') 2000 FORMAT (F9.4) C 300 IF (IC.NE. 1) GO TO 500 WRITE (06,1000) READ (05,2000) RA(IV) WRITE (06,111) IV,RA(IV) 111 FORMAT (' — VARS ARRAY —-',I6,F10.4) GO TO 900 C 500 IF (IC.NE. 2) GO TO 800 WRITE (06,1000) IF (IV.GT. 7) GO TO 700 READ (05,2000) RB1(IV) WRITE (06,222) IV,RBI(IV) 222 FORMAT (' — PAR ARRAY —-',IG,F10.4) GO TO 900 C 700 IN = IV-7 READ (05,2000) RB2(IN) WRITE (06,222) IN,RB2(IN) C 800 WRITE (06,1002) C 900 CONTINUS RETURN END