Some Studies on Modeling and Simulation of a Disc Brake System for
Squeal Prediction
Except where reference is made to the work of others, the work described in this thesis
is my own or was done in collaboration with my advisory committee. This thesis does
not include proprietary or classified information.
Nitesh Ahuja
Certificate of Approval:
David G. Beale
Professor
Mechanical Engineering
P.K.Raju, Chair
Thomas Walter Professor
Mechanical Engineering
Dan B. Marghitu
Professor
Mechanical Engineering
Stephen L. McFarland
Acting Dean
Graduate School
Some Studies on Modeling and Simulation of a Disc Brake System for
Squeal Prediction
Nitesh Ahuja
A Thesis
Submitted to
the Graduate Faculty of
Auburn University
in Partial Fulfillment of the
Requirements for the
Degree of
Master of Science
Auburn, Alabama
August 7, 2005
Some Studies on Modeling and Simulation of a Disc Brake System for
Squeal Prediction
Nitesh Ahuja
Permission is granted to Auburn University to make copies of this thesis at its
discretion, upon the request of individuals or institutions and at their
expense. The author reserves all publication rights.
Signature of Author
Date of Graduation
iii
Vita
Nitesh Ahuja was born in Faizabad, Uttar Pradesh, India on January 12, 1978.
He graduated with a Bachelor of Engineering degree in Mechanical Engineering from
Dayananda Sagar College of Engineering, Bangalore University, Bangalore, India in De-
cember 2001. He joined the Masters program in the Department of Mechanical Engi-
neering at Auburn University in August 2003.
iv
Thesis Abstract
Some Studies on Modeling and Simulation of a Disc Brake System for
Squeal Prediction
Nitesh Ahuja
Master of Science, August 7, 2005
(B.E., Bangalore University, 2002)
92 Typed Pages
Directed by P.K.Raju
Most mechanical systems that involve a sliding contact experience some degree of
friction induced self-excited vibration. A well known example is brake squeal which is
a top concern for automotive companies, who must spend huge amounts of money on
warranty repairs. Brake squeal occurs above 1k Hz and below 20k Hz. The phenomena
that cause brake squeal is not thoroughly understood, although a great deal of research
has been dedicated to developing a better understanding of the complex problem of brake
dynamics. This study developed a displacement dynamic equation of motion for brake
rotor?s, assuming the rotor?s to be annular plates with the same modal characteristics.
An analytical approach was used to solve the equation of motion of these annular discs
for free vibration. The Eigen values obtained through the analytical approach were then
used to obtain the natural frequencies. Experimental modal testing was also performed
to determine the natural frequencies of disc and pad with free boundary conditions.
The natural frequencies of the annular disc were also obtained using the finite dif-
ference method. An algorithm was developed to describe the dynamic response of the
v
annular disc for both a constant forcing function and a time dependent forcing function
using the central difference method. The stability of the algorithm was also examined
in order to obtain the time step. The importance of the finite difference technique lies
in its ability to visualize the results and provide deeper insights into the dynamics of
brake rotors. Computer software in C programming and Matlab was developed to solve
the differential equation with the forcing function. Another advantage of this method
is that any type of forcing function which is practically possible can be applied and its
response calculated. The results revealed that an annular disc behaves non-linearly.
Finite element analysis was then performed to determine the dynamic character-
istics (natural frequency and mode shapes) and the response of the brake model. The
advantage of using this technique is that frictional forces are also taken into account. The
contact analysis applied a special type of element called CONTAC49 for 3-D surfaces.
In this thesis a four degree of freedom mathematical model of the brake system
based on friction characteristics and geometric non-linearity capable of representing the
brake dynamics is presented. Friction induced vibration was studied by using a four
degree of freedom model consisting of a disc and a pad moving in both the in-plane and
transverse directions. This system takes the form of a mathematical model with two
coupled equations and two uncoupled equation. The set of four differential equations
can then be solved using Runge Kutta?s JBM method.
vi
Acknowledgments
The author would like to thank his advisor Dr.P.K.Raju,the Thomas Walter Pro-
fessor and Chair, Mechanical Engineering Department for his guidance and support
throughout this research and especially during the completion of this thesis.
The author would also like to thank Dr.David G. Beale, Professor, Mechanical Engineer-
ing Department and Dr. Dan B. Marghitu, Associate Professor, Mechanical Engineering
Department, for serving as committee members. Finally, the author thanks his parents
and friends for their support and help.
vii
Style manual or journal used LATEX: A Document Preparation System by Leslie
Lamport (together with the style known as ?aums?).
Computer software used The document preparation package TEX (specifically
LATEX) together with the departmental style-file aums.sty. The images and plots were
generated using MATLAB 7.0 and Microsoft Office Excel 2003.
viii
Table of Contents
List of Figures xi
List of Tables xiii
1 Introduction 1
1.1 Scope of the Thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 An Overview of the Thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
2 Background 4
2.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.2 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
3 Mathematical Modeling of A Disc Brake 10
3.1 Free Vibration of Circular Disc Without Damping : An Analytical Approach 10
3.1.1 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.2 Experimental Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.3 Dynamic Equation of Flexible Annular Disc . . . . . . . . . . . . . . . . . 18
3.4 Finite Difference Discretization . . . . . . . . . . . . . . . . . . . . . . . . 22
3.5 Initial and Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . 26
3.6 Stability of the Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.7 Extraction of Eigenvalues and Eigenvectors: Finite Element Approach . . 35
3.7.1 Block Lanczos Method . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.7.2 Contact Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
4 Theoretical Modeling of Disc-Brake System 42
4.1 Lagrange Formulation and Equation of Motion . . . . . . . . . . . . . . . 42
4.1.1 Friction Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
4.2 Numerical Simulation and Stability Analysis . . . . . . . . . . . . . . . . . 46
5 Results and Conclusions 50
5.1 Summary of Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
5.1.1 Analytical Method . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
5.2 Finite Difference Method . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
5.2.1 Case I . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
5.2.2 Case II . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
5.3 Finite Element Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
5.4 Four Degrees of Freedom Model . . . . . . . . . . . . . . . . . . . . . . . . 54
5.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
ix
Bibliography 57
Appendices 59
A C-program to solve nonlinear differential equation using finite dif-
ference method 60
B Matlab program to predict the response of the disc 71
C C-program to solve Coupled Differential Equation 72
x
List of Figures
2.1 Disc Brake System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
3.1 Annular disc brake . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
3.2 Figure shows SC-2345, accelerometers, impact hammer, foam support,
PCMCIA card and computer system . . . . . . . . . . . . . . . . . . . . . 19
3.3 Disc modal peaks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.4 Pad modal peaks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.5 Representation of grid points using a 13 node mesh in FDM . . . . . . . . 23
3.6 Response due to a force applied at node 4 of 5X5 grid points . . . . . . . 28
3.7 Response due to a force applied at node 16 of 5X5 grid points . . . . . . . 28
3.8 Response due to a force applied at node 22 of 5X5 grid points . . . . . . . 29
3.9 Response due to a force applied at nodes 2,7,12 and 22 of 5X5 grid points 29
3.10 Response due to a force applied at nodes 3,8,13,18 and 23 of 5X5 grid points 30
3.11 Response due to a force applied at nodes 18,19,23 and 24 of 5X5 grid points 30
3.12 Forcing Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.13 Dynamic response in transverse direction due to force applied at point 20
of 10X5 grid points . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.14 Dynamic response in transverse direction due to force applied at point 23
of 10X5 grid points . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.15 Dynamic response in transverse direction due to force applied at point 28
of 10X5 grid points . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.16 Another view of the meshed brake rotor model . . . . . . . . . . . . . . . 39
3.17 Equivalent disc model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
xi
3.18 Modeshape of a disc at natural frequency of 1205 Hz . . . . . . . . . . . . 40
3.19 Modeshape of a disc at natural frequency of 1840 Hz . . . . . . . . . . . . 40
3.20 Modeshape of a disc at natural frequency of 2937 Hz . . . . . . . . . . . . 40
3.21 The generation of contact element Contac49 at the disc and pad interface
in ANSYS 8.0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.22 Pressurized brake disc model in ANSYS 8.0 . . . . . . . . . . . . . . . . . 41
4.1 Four degrees of freedom model . . . . . . . . . . . . . . . . . . . . . . . . 44
4.2 (a)Motion & (b)Time history of the disc (Force N=1 and b=0.1)respectively 48
4.3 (a)Motion & (b)Time history of the disc (Force N=1 and b=1)respectively 49
4.4 (a)Motion & (b)Time history of the disc (Force N=1 and
us=0.01)respectively . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
4.5 (a)Motion & (b)Time history of the disc respectively (Force N=1 and
us=0.1) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
5.1 Response due to a force applied at node 17 of 5X5 grid points . . . . . . . 52
5.2 Forcing Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
5.3 Dynamic response in transverse direction due to a force applied at point
17 of 10X5 grid points . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
5.4 Nodal displacement of the disc in the z direction . . . . . . . . . . . . . . 54
5.5 Reponse in the z direction with repect to time . . . . . . . . . . . . . . . . 55
xii
List of Tables
3.1 Parameters of the annular plate . . . . . . . . . . . . . . . . . . . . . . . . 16
3.2 Natural frequencies of the disc obtained by analytical approach . . . . . . 17
3.3 Natural frequencies (Hz) of the disc found experimentally . . . . . . . . . 18
3.4 Natural frequencies of the disc obtained by the finite difference method . 27
3.5 Natural frequencies (Hz) of the disc using finite element analysis . . . . . 39
5.1 Natural frequencies of the disc obtained through all methods . . . . . . . 50
5.2 Percentage deviation of natural frequencies of the disc with respect to the
experimental method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
xiii
Chapter 1
Introduction
Brakes are one of the most vital safety features in motorized vehicles. However, Au-
tomotive brake noise has become an increasingly important issue for automobile manu-
facturers. Today, most manufacturers attempt to reduce the brake noise by changing the
geometry, carefully selecting the friction materials and mounting vibration shims at the
backplate of the pads. However, the mechanism that causes this noise is not yet properly
understood. This report presents a new analysis methodology which has the potential
to improve the efficiency and accuracy of the design process for disc-brake systems. This
process will allow a more accurate determination of the variable parameters influence
the brake noise and that reduces the time needed to complete the design.
1.1 Scope of the Thesis
Structural dynamics and friction play an important role in disc-brake dynamics.
For the brake system studied in this thesis, the rotor was the major cause of squeal
and radiator of sound. The natural frequency of a rotor can be obtained analytically by
solving a differential equation with boundary conditions applied. The aim of the work
presented in this thesis has been to determine natural frequency of the disc theoretically.
The thickness of the annular disc is assumed to act in such a way that its natural fre-
quencies match the actual natural frequency of the rotor of the vehicle. After verifying
the theoretical model, the dynamics of a disc subjected to nodal forces were studied. The
natural frequencies obtained theoretically were verified by comparison with the natural
1
frequency obtained experimentally. Various type of static loading and dynamic loading
were applied at the different nodes and the response of the disc predicted. Finite Ele-
ment Analysis was also performed to understand the behavior as disc and brake pads
interact. To test the behavior of the model, a contact element was generated between
the two surfaces using the three dimensional contact element contac49 in ANSYS8.0.
A four degree of freedom mathematical model was also developed to simulate a disc
brake system with a cubic non-linearity. Using this model, the interaction between the
pad and rotor was investigated. Non-linear analysis was also performed to demonstrate
the unstable spiral and limit cycles in the phase space. The results showed that one of
the instability conditions depend on the initial velocity of the disc. It also provided new
insights that will enable us to understand contact and friction phenomena in the brake
pad and disk coupling.
1.2 An Overview of the Thesis
This thesis presents an in-depth study of the dynamic behavior of a brake disc using
finite difference method.
A review of the literature on disc-brake noise is presented in Chapter 2. The natural
frequencies of an annular plate are found analytically, which match the actual brake rotor
and a mathematical equation is developed for the annular plate. Due to the complex-
ity of the mathematical model, it must be solved using the Finite Difference Method.
Chapter 4 presents the theoretical modeling of a disc brake system cosidering friction and
describes the experimental work performed in order to obtain a better understanding of
2
the squealing phenomena. Chapter 5 discusses the results and presents the conclusions
that can be drawn as a result of the analysis.
3
Chapter 2
Background
2.1 Overview
Although, Automotive brakes have improved significantly during the last decade,
the engine power has also increased, so the energy to be dissipated in a brake is several
times larger than the power developed by the engine. However, not only the power but
also the expectations of comfort by drivers and passengers have also greatly increased
over the last decade. This means that noise levels, and in particular brake noise, which
were acceptable 20 or 30 years ago are no longer acceptable. Brake squeal is now a ma-
jor problem in vehicles leading to substantial losses for automobile companies. The disc
brake was invented in 1902 by Frederick William Lanchester in England[1]. Nowadays,
the vibration and noise produced in disc-brake systems cause discomfort for passengers
and raise their concerns over the reliablity of the vehicle. Akay[2] recently reviewed
friction-induced noise, including disk-brake squeal. He quoted an industrial source who
gave an estimate of the warranty costs due to noise, harshness and vibration (together
known as NVH problems), including disk-brake squeal, as costing the automotive indus-
try U.S. $ 1 billion in North America alone.
Brake noise is usually classified in three categories : low frequency noise, low fre-
quency squeal and high frequency squeal. Low frequency brake noise typically occurs
between 100 Hz and 1 KHz and is referred to as grunt, groan, grind and moan. It
originates in the excitation of the friction material at the rotor and lining interface. In
4
Figure 2.1: Disc Brake System
disk brakes, squeal can occur when the brake pads contact the rotor while the vehicle is
moving at low speeds, causing a vibration that manifests as an annoying high-pitched
squeal. Low frequency squeal occurs in the frequency range above 1 kHz and below
the first circumferential mode of the rotor[3]. This type of noise is caused frictional
excitation and the coupling of more than two modes of various component of the brake
system. This coupling produces optimum conditions for brake squeal. High frequency
squeal is a particularly troublesome problem in brake development and typically occurs
at frequencies above 5 kHz. It has been experimentally demonstrated that the frequency
of squeal coincides with the circumferential resonance fequencies of the rotor disc[3].
A car disc brake system consists of a rotating disc, non-rotating friction pads, carrier
bracket, calliper, yoke, piston,sealing ring, guide pins and dustboot. The pads are loosely
housed in the calliper and located by the carrier bracket. The calliper itself is allowed
5
to slide fairly freelyalong the two mounting guide pins in a floating calliper design. The
disc is bolted to the car wheel and thus rotates at the same speed as the wheel. When
the disk brake is applied, the two pads are brought in contact with the rotating disc sur-
faces. Most of the kinetic energy of the vehicle is converted to heat through friction but
a small part of it is converted into sound energy and generates noise.[4]. Optimizing this
automotive noise behavior is an important aspect in new vehicle designs. The purpose
of a brake is to attenuate effectively and evenly the kinetic energy of a vehicle. This is
done by converting it into heat energy which is generated due to the friction between
the brake disk and the pads. Under certain operating parameters, during the braking
process, high frequency vibrations of the brake disc can be excited. These significant
structural vibrations of the brake, which have no influence on the actual braking effect,
are heard as an undesired, acoustically perceptible squealing of brakes. The out-of-plane
vibrations are physically responsible for the noise emitted from the brake disk when
squealing occurs. In this case, the surface of the friction ring of the brake disc vibrates
at right angles to the disc plane. The decelerating braking power between pads and
brake disc has an effect in the disc plane, exciting in-plane vibrations.[5]
2.2 Literature Review
North [7] observed that the first concerted effort to study car disk brake squeal was
made at the Motor Industrial Research Association (MIRA) in the U.K. in the 1950s.
In 1970, a study on the vibration of a disc subject to a moving load was initiated. Mote
6
investigated the vibration of a disc under a pointwise load and showed that the instabil-
ity is likely to occur in the supercritical range [6].
Disc brake squeal occurs when a system experiences a very large amplitude vibra-
tion. There are two theories that attempt to explain why this occurs. The first theory
states that a negative friction-velocity gradient leads to stick-slip vibration, which has
been identified as one of the major causes of brake squeal, and it occurs due to non-
linearity of the frictional force-velocity characteristics. The second theory believes that
the high level of vibration results from geometric instabilities in the brake assembly.
Both theories attribute the brake system vibration and the accompanying audible noise
to variable friction forces at the lining-rotor interface [12]. Yuan classified these mech-
anisms into two groups [13]. The friction coefficient tends to decrease as sliding starts.
This is the basis of the negative slope of friction versus relative speed, which was first
put forward as a possible mechanism responsible for generating brake squeal by Mills
and was later supported by experimental results [14]. However, North [7] thought that
it(alone) was not likely to account for squeal. The difference between the static friction
coefficient and the kinetic friction coefficient may cause stick-slip vibration, which was
first introduced to disk brake squeal studies by Mills [14]. North used a two-degree-of-
freedom model of a car disc brake to show that the different friction forces on either side
of the disc brake could lead to instability. In another study, a slider experiencing stick-
slip vibration on a vibrating disk was shown to display a rich dynamic behavior [15].
Ouyang, Cao, Mottershead and Treyde [16] believe that a good disc brake model should
include a plausible squeal mechanism, a proper treatment of the frictional contact and an
accurate description of the system dynamics. The disc is presented as an annular plate
7
and the stationary components (pad, calliper, carrier and mounting pins) are described
using a finite element model, this constituting a combined numerical-analytical model.
The fricition forces acting on the surface of the disk generate pointwise bending couples
that move around the disc. Thus, the whole disk brake system must be treated as a
proper moving load problem through the moving contact at the disk/pads interface with
friction [16]. Ouyang and Mottershead [17] investigate the instability of the transverse
vibration of a disc excited by two co-rotating sliders on both sides of the disc, consid-
ering each slider as a mass-spring-damper system traveling at the same constant speed
around the disc. In this model, Frictional forces acting at the contact interface between
the disc and each of the two sliders. The bending couple and normal forces produced by
the rotating sliders are moving loads and are found to result in dynamic instability [17].
Chowdhary, Bajaj and Krousgrill [18] predict one of the mechanisms behind the disc
brake model, where the brake rotor is represented by a thin plate of equivalent modal
characterisitics and the backing plate are modeled as thin annular section plates using
the Rayleigh-Ritz approach. The two structural models are then coupled using linear
elastic springs and Coulomb friction at the interface, and langrangian approach is used
to derive the equations of motion of the coupled system [18].
Whatever the modeling method adopted to solve this problem, linear analysis is
commonly used to determine the stability of the system for various parametric condi-
tions. By performing stability analysis, the real parts of the eigenvalues can be extracted
to construct the characteristic equation of the system. If the real parts are positive then
the system is unstable, and if they are negative then the system is stable. In a disc brake
system, friction force is the main cause for the excitation of vibration and is unevenly
8
distributed force between the brake disc and the brake pad on either sides of the rotor.
Research on the vibrational instabilities which lead to brake squeal has been con-
ducted for more than fifty years. Today?s better understanding of the vibrational insta-
bilities that produce brake squeal has required a great deal of detailed and sophisticated
analysis and experimental research due to the complexity of a brake system. The ability
to design a brake system and identify its causes of squealing in advance would offer a
very helpful design tool. The literature also suggests that different brake systems radiate
sound in very different ways. Murakami(1984) studied a brake system and found that
the brake calliper and pad played major role in brake squeal. Five years later, Nishiwaki
conducted a research on a different brake system and concluded that rotor was respon-
sible for the most noise. McDaniel, Moore, Chen and Clarke [19] presented a study of
acoustic radiation from the stationary brake system. The objective of this experiment
was to better understand the acoustic radiation from squealing brake systems. The re-
searchers found that the great majority of squeal mechanisms occur due to the resonant
behavior of the operating brake system. They also presented an analysis that equated
the natural frequencies and modes of mechanically-excited stationary brake systems to
those of an operating brake system. Acoustic radiation efficiencies and intensities of the
modes were computed by importing experimentally measured velocities into a BEM soft-
ware package, which revealed that for a particular brake system, the highest radiation
efficiency occurred at frequencies above 2-3kHz [19].
9
Chapter 3
Mathematical Modeling of A Disc Brake
3.1 Free Vibration of Circular Disc Without Damping : An Analytical Ap-
proach
It is difficult to write a differential equation for a disc-brake model, which has a
very complicated geometry. In order to develop a mathematical model for a disc-brake
system, it is first necessary to reduce the disc into an annular plate with the same natural
frequency and mode shapes. In this case, the first nine natural frequencies are matched
with the actual rotor, which is achieved by changing the thickness of the annular plate.
The displacement equation of motion for circular plate without damping can be
written as:
D?4w +?h?w = p(r,?,t) (3.1)
where
?2 = ?
2
?r2 +
1
r
?
?r +
1
r2
?2
??2
D = E ?h
3t
12?(1??2)
10
Figure 3.1: Annular disc brake
E is the Elastic modulus, h is the thickness,? is the mass desity and ? is the Poisson?s
ratio. Considering the case of free vibration when the forcing function at the right hand
side i.e. p(r,?,0)=0 and assuming the boundary conditions to be homogeneous, the above
equation can be written as
D?4w +?h?w = 0 (3.2)
The homogeneous boundary condition can take one of the forms
Vn = 0 (3.3)
or ,
Mnn = 0 (3.4)
11
Thus, either
w = 0 (3.5)
or,
?w/?n = 0 (3.6)
Let us assume that the free vibrations are characterized by
w(r,?,t) = W(i)(r,?)cos(?it),i = 1,2,3... (3.7)
w(r,?,t) = W(j)(r,?)cos(?jt),j = 1,2,3... (3.8)
where ?i is the natural frequency and Wi(r,?) is its associated mode shape, then upon
substitution of equation 3.7 and 3.8 into 3.2 [20], we obtain
D?4W(i) = ?h?2iW(i) (3.9)
D?4W(j) = ?h?2jW(j) (3.10)
The equations 3.9, 3.10 resemble static equilibrium equations and therefore the
Betti-Rayleigh reciprocal theorem [20] is applicable if it is set as pi=? h ?2i Wi and pj=?
h ?2j Wj. Assuming that the present analysis is restricted to homogeneous boundary
conditions, the application of the Betti-Rayleigh reciprocal theorem to the result obtained
12
in 3.10 and 3.11, gives
integraldisplay
A
?h?2iW(i)W(j)dA =
integraldisplay
A
?h?2jW(j)W(i)dA
(?2i ??2j)
integraldisplay
A
?h?2jW(i)W(j)dA = 0 (3.11)
If ?2i negationslash= ?2j , then
integraldisplay
A
?h?2jW(i)W(j)dA = 0 (3.12)
3.1.1 Boundary Conditions
Consider the axisymmetric, harmonic vibration of the disc clamped at radius r1.
If the plate deflection is referred to in polar coordinates. We can assume the axial
symmetry, w=w(r), where r is the radial coordinate. The mode shapes W(i) must satisfy
the boundary condition, as the disc is fixed or clamped at the inner boundary.
At the inner boundary i.e at r=r1 we have,
W(i)(r) = ?W
(i)
?r = 0 (3.13)
Therefore,
(?2 +?2i)(?2 ??2i)W(i) = 0 (3.14)
where
?2i = ?i(?hD )1/2 (3.15)
13
The boundary condition at the outer edge of the disc i.e at r=r2 is assumed to be
free,
Therefore, the bending moment and shear force at the outer edge is zero.
Mathematically, this can be written as follows:
Mr |(r=r2)= ?D[?
2W
?r2 +?(
1
r
?W
?r +
1
r2
?2W
??2 )] = 0 (3.16)
Vr |(r=r2)= ?[D ??r(?2W)+(1??) ?r??(1r ?
2W
?r?? ?
1
r2
?W
?? )] = 0 (3.17)
Since ?2i > 0,the solution of equation 3.14 can be written in the form,
W(r,?) = (C1cos(?)+C2sin(?))(AiBesselJ(0,?ir)+BiBesselY(0,?ir)
+CiBesselI(0,?ir)+DiBesselK(0,?ir)) (3.18)
where J0 and K0 are Bessel functions, and I0 and K0 are modified Bessel functions, all of
order zero. The above assumed solution should also satisfy all the boundary conditions
as follows:
W(r=r1) = AiBesselJ(0,?ir1)+BiBesselY(0,?ir1)+CiBesselI(0,?ir1)
+DiBesselK(0,?ir1) = 0 (3.19)
?W
?r (r=r1) = ?AiBesselJ(0,?ir1)?i ?BiBesselY(0,?ir1)?i +CiBesselI(0,?ir1)?i
?DiBesselK(0,?ir1)?i = 0
14
(3.20)
The bending moment at the outer edge is given by the following equation:
M(r=r2) = ?Ai[BesselJ(0,?ir2)? BesselJ(1,?ir2)?
ir2
]?2i ?Bi[BesselY(0,?ir2)
?BesselY(1,?ir2)?
ir2
]?2i +Ci[BesselI(0,?ir2)? BesselI(1,?ir2)?
ir2
]?2i
?Di[BesselK(0,?ir2)? BesselK(1,?ir2)?
ir2
]?2i + 1r
2
[?(?AiBessel(1,?ir2)?i
?BiBesselY(1,?ir2)?i +CiBesselI(1,?ir2)?i ?DiBesselK(1,?ir2)?i)]
+ ?r2[?AiBesselJ(0,?ir2)?BiBesselY(0,?ir2)?CiBesselI(0,?ir2)
?DiBesselK(0,?ir2)] = 0
(3.21)
V(r=r2) = (1??)?ir2 [AiBesselJ(1,?ir)+BiBesselJ(1,?ir)?CiBesselJ(1,?ir)
+DiBesselJ(1,?ir)]+ (1??)r3 [AiBesselJ(0,?ir)+BiBesselJ(0,?ir)
?CiBesselJ(0,?ir)+DiBesselJ(0,?ir)]?Ai?2i[?BesselJ(1,?ir)?i
+2BesselJ(1,?ir)?
ir2
? BesselJ(0,?ir)r ]?Bi?2i[?BesselY(1,?ir)?i
+2BesselY(1,?ir)?
ir2
? BesselY(0,?ir)r ]+Ci?2i[BesselI(1,?ir)?i
+2BesselI(1,?ir)?
ir2
? BesselI(0,?ir)r ]?Di?2i[BesselK(1,?ir)?i
+2BesselK(1,?ir)?
ir2
+ BesselK(0,?ir)r ]?Ai?2i[BesselJ(0,?ir)?ir?BesselJ(1,?ir)?
ir2
]
?Bi?2i[BesselY(0,?ir)?ir?BesselY(1,?ir)?
ir2
]+Ci?2i[BesselI(0,?ir)?ir?BesselI(1,?ir)?
ir2
]
?Di?2i[?BesselK(0,?ir)?ir?BesselK(1,?ir)?
ir2
]+AiBesselJ(1,?ir)?ir2
15
Table 3.1: Parameters of the annular plate
Parameter Value
Thickness of disc,ht 0.012 m
Inner radius of disc, r1 0.045 m
Outer radius of disc, r2 0.133 m
Young?s Modulus, E 1.2E11 Pa
Mass Density, ? 7200 Kg/m3
Poisson?s ratio, ? 0.211
+BiBesselY(1,?ir)?ir2 ?CiBesselI(1,?ir)?ir2 +DiBesselK(1,?ir)?ir2
?AiBesselJ(0,?ir)r2 ?BiBesselY(0,?ir)r2 ?CiBesselI(0,?ir)r2 ?DiBesselK(0,?ir)r2 = 0
(3.22)
Solving equations 3.19, 3.20, 3.21 and 3.22 for Ai, Bi, Ci and Di and substituting those
values of Ai, Bi, Ci and Di into equation 3.19, and again solving using the Maple 10
software for ? at r=r1. From equation 3.14 , the natural frequencies of vibration are
obtained by finding the roots of the equation 3.19 as follows:
?2i = ?i(?hD )1/2
?i = ?
2i
(?hD )1/2 (3.23)
The parameters given in table 3.1 were used to calculate the natural frequencies of
the disc. The natural frequencies obtained analytically listed in table 3.2
16
Table 3.2: Natural frequencies of the disc obtained by analytical approach
Serial Number Value of ? Natural Frequency (Hz)
1 7.54 979
2 8.12 1135
3 8.39 1212
4 10.3 1827
5 12.98 2902
6 15.831 4316
7 19.1 6282
8 20.27 7075
9 20.45 7201
3.2 Experimental Setup
The pad and the disc were set into vibration at a fixed-point by an ICP Impulse
Force Test Hammer model 086D80, provided by PCB Piezotronics. The responses were
measured at a number of points on the surface of the structure using a piezoelectric ICP
accelerometer model number 352C15 provided by PCB Piezotronics. Two SCC-ACC01
accelerometer input modules with a fixed gain of 2 provided a 4mA current source for
the accelerometer and impulse hammer. The 2 modules accepted input signals from the
accelerometer and impulse hammer, amplified them and coupled these signals to the
data acquisition system via an SC-2345 Signal Conditioning with configurable connec-
tors. The SC-2345 carried signals from 68-pin E Series acquisition (DAQ) devices. It
was used with SCC Series modules and a shielded 68-pin cable and offered easy-to-use,
low-noise signal conditioning on a per-channel basis.
In this experiment, the signal acquisition was performed using the NI DAQ Card-6036E,
16 inputs/2 outputs, 200 kS/s and LabView software for A/D data conversion, data
processing, analysis and archival.Instruments for data recording stored data in an exist-
ing file. The diagram also contained instruments for windowing, spectrum units selection
17
Table 3.3: Natural frequencies (Hz) of the disc found experimentally
Serial Number Natural Frequency(Hz)
1 998
2 1161
3 1205
4 1840
5 2937
6 4380
7 6212
8 7026
9 7295
and immediate display of time waveforms and frequency spectra.
After the data was properly acquired and analyzed, ME?scope VES 4.0 software, Vi-
sual Modal Pro package for experimental modal testing was applied. ME?scope VES 4.0
was developed by Vibrant Technology and consists of a set of post-test analysis software
that can be used to analyze and document the dynamic behavior of mechanical struc-
tures and machines. The Visual Modal Pro package imports the FRF measurements and
generates estimates by curve fitting the modal parameters, namely the natural frequen-
cies, damping and mode shapes.
3.3 Dynamic Equation of Flexible Annular Disc
Spinning annular discs are commonly found in mechanical as well as in electrical
engineering applications such as turbine rotors, computer disc memory units, vehicle
disc-brake etc. Therefore, their dynamic analysis is of major interest to researchers [21].
It is essential to understand the flexible nature of an annular disc in order to construct
18
Figure 3.2: Figure shows SC-2345, accelerometers, impact hammer, foam support, PCM-
CIA card and computer system
Figure 3.3: Disc modal peaks
19
Figure 3.4: Pad modal peaks
an accurate mathematical model for the disc-brake system, which would be helpful in
controling the vibration of the disc efficiently. This kind of model can be developed
using partial differential equations to describe the dynamic behavior of a flexible disc.
The partial differential equation can then be solved by discretizing it into a finite set of
ordinary differential equations. Finite element analysis has also been used to understand
the dynamic nature of the plate. For example, A new symmetric and positive definite
boundary element formulation for lateral vibration of plates was introduced by Dav and
Milazzo [22].
Consider an annular plate that is spinning at an angular speed of ? and is transverly
in contact with a stationary oscillating unit at the point (rp,0), where (r,?) is a polar
20
coordinate system fixed in space. The plate is clamped or fixed at the inner edge and
free at the outer boundary.
D?4w +c(?w?t )+?h(?
2w
?t2 ) = F(r,?,t)+ht[
1
r
?
?r(r?rr
?w
?r )+
1
r2???
?2w
??2 ] (3.24)
where :
?4 = Biharmonic operator,
w = lateral deflection in the z direction, a function of r, ? and t.
? = mass density per unit area,
D = Eh3t(1??2) flexural rigidity,
ht = thickness of the plate,
? = Poisson?s ratio,
E = Young?s Modulus,
c = viscous damping,
?2w
?t2 = acceleration in the z direction,
?rr and ??? = Initial in-plane stresses induced by rotation.
The above equation represents the equation of motion characterizing the behavior of a
flexible annular disc in lateral motion. The solution of this equation can be obtained by
considering the corresponding boundary conditions. Although the equation contains a
term for damping, this will be neglected damping at this stage when solving the above
equation. Also,the initial stresses induced due to rotation can be neglected as their effect
will be negligible compared to the other stresses present. The above equation can thus
21
be rewritten as:
D?4w +?h(?2w/?t2) = F(r,?,t) (3.25)
A numerical method approach based on finite differences can be applied to simulate the
characterising behavior of an annular disc. This above equation can also be written in
expanded form as follows:
D[?
4w
?r4 +
2
r
?3w
?r3 +
1
r3
?w
?r +
2
r3
?3w
??2?r +
2
r2
?4w
??2?2r
+ 4r4 ?
2w
??2 +
1
r2
?2w
?r2 +
1
r4
?4w
??4 ]+?h
?2w
?t2 = 0 (3.26)
3.4 Finite Difference Discretization
The difference between the analytical solution and finite difference method is that
the solution is obtained only at discrete points. Therefore, dividing the region into
small regions and assigning each region a reference index will be provide very accurate
results. Solving the circular boundary problem solved using the finite difference approach
defined in polar co-ordinates (r,?) is more convenient than using cartesian co-ordinates as
it avoids the use of awkward differentiation formulae near the curved boundary. Defining
the mesh in the r-? plane by the points of intersection of the circles r=r1 + i*h, where
i=0, 1, 2, 3 ..., r1 is the inner radius of the disc and h=(r2-r1)/M, M is the number of
circles defined for the mesh and the straight lines ?= j * l where j=0, 1, 2, 3.... and
l=(?2-?1)/N , where N is the number of divisions. These lines and circles are called grid
lines and their intersections are referred to as mesh points of the grid or nodal points.
The distance between two nodal points is the grid size. In some cases, the accuracy of
22
Figure 3.5: Representation of grid points using a 13 node mesh in FDM
the solution can be improved by reducing the mesh size. Figure 3.5 shows a nodal point
at the center of a small shaded region that has the index (i, j) and its neighbouring points
have node indices (i+1, j), (i-1, j),(i, j+1), (i, j-1), (i+1, j+1) and so on. In the case of
a flexible annular plate, a three dimensional coordinate system is considered. Here, the
third dimension is time, t. The time axis is represented with reference index k, where
t=k*?t. For each nodal point at the interior of the grid (ri, ?j, tk) where i=0, 1, 2, ...m,
j=0, 1, 2, ...n, k=1, 2, ....q. A Taylor series expansion is used to generate the central
finite difference for the partial derivative terms of the lateral deflection, w(r, ?, t)= wi,j,k
The above equation for the free vibration of the disc without damping can be ex-
panded using the central difference method and can be written as,
[?
4w
?r4 ]i,j,k =
1
(r1+i?h)4[wi+2,j,k ?4wi+1,j,k +6wi,j,k ?4wi?1,j,k +wi?2,j,k] (3.27)
1
r4[
?4w
??4 ]i,j,k =
1
(l)4(r1+i?h)4[wi,j+2,k?4wi,j+1,k+6wi,j,k?4wi,j?1,k+wi,j?2,k] (3.28)
23
2
r[
?3w
?r4 ]i,j,k =
2
h3 ?(r1+i?h)[wi+2,j,k ?2wi+1,j,k +2wi?1,j,k ?wi?2,j,k] (3.29)
1
r3[
?w
?r ]i,j,k =
1
h?(r1+i?h)3[wi+1,j,k ?wi?1,j,k] (3.30)
2
r3[
?3w
??2?r]i,j,k =
2
(l)2 ?h?(r1+i?h)3[wi+1,j+1,k ?wi?1,j+1,k ?2wi+1,j,k
+2wi?1,j,k +wi+1,j?1,k ?wi?1,j?1,k] (3.31)
2
r2[
?4w
??2?r2]i,j,k =
2
(l)2 ?h2 ?(r1+i?h)2[4wi,j,k ?2wi+1,j,k
?2wi?1,j,k ?2wi,j+1,k ?2wi,j?1,k +wi+1,j+1,k
+wi+1,j?1,k +wi+1,j+1,k +wi+1,j+1,k] (3.32)
4
r4[
?2w
??2 ]i,j,k =
4
(l)2 ?(r1+i?h)4[wi,j+1,k ?2wi,j,k +wi,j?1,k] (3.33)
1
r2[
?2w
?r2 ]i,j,k =
1
(h)2 ?(r1+i?h)2[wi+1,j,k ?2wi,j,k +wi?1,j,k] (3.34)
?h[?
2w
?t2 ]i,j,k =
?h
?t2[wi,j,k+1 ?2wi,j,k +wi,j,k?1] (3.35)
Substituting ?4w?r4 , 1r4 ?4w??4 , 2r[?3w?r4 ], 1r3[?w?r ], 2r3[ ?3w??2?r], 2r2[ ?4w??2?r2], 4r4[?2w??2 ], 1r2[?2w?r2 ],
?h[?2w?t2 ] from equations (3.25-3.33) into equation (2.16) and simplifying yields:
wi,j,k+1 = ?D??t
2
??ht (Pwi,j,k +Qwi+1,j,k +Rwi?1,j,k +S(wi,j+1,k +wi,j?1,k)
+T(wi+1,j+1,k +wi+1,j?1,k)+U(wi?1,j+1,k +wi?1,j?1,k)
+V(wi+2,j,k +wi?2,j,k)+W(wi,j+2,k +wi,j?2,k))+2wi,j,k
?wi,j,k?1 + ?t
2F(r,?,t)
?ht (3.36)
24
where
P = [ 6h4 + 6((r1+i?h)4 ?(??4) + 8(r1+i?h)2 ?(h2)?(??2)
? 8(r1+i?h)4 ?(??2) ? 2(r1+i?h)2 ?(h2)],
Q = [? 4h4 ? 4(r1+i?h)?(h)3 + 1(r1+i?h)3 ?h ? 4r1+i?h)3 ?h?(??2)
? 4(r1+i?h)2 ?(h2)?(??2) + 1(r1+i?h)2 ?(h2)]? 8(r1+i?h)4 ?(??2)
? 2(r1+i?h)2 ?(h2)],
R = [? 4h4 + 4(r1+i?h)?(h)3 ? 1(r1+i?h)3 ?h + 4r1+i?h)3 ?h?(??2)
? 4(r1+i?h)2 ?(h2)?(??2) + 1(r1+i?h)2 ?(h2)]? 8(r1+i?h)4 ?(??2)
? 2(r1+i?h)2 ?(h2)],
S = [? 4((r1+i?h)4 ?(??4) ? 4(r1+i?h)2 ?(h2)?(??2) + 4(r1+i?h)4 ?(??2)]
? 8(r1+i?h)4 ?(??2) ? 2(r1+i?h)2 ?(h2)],
T = [ 2((r1+i?h)3 ?h?(??2) + 2(r1+i?h)2 ?(h2)?(??2)]? 8(r1+i?h)4 ?(??2)
? 2(r1+i?h)2 ?(h2)],
25
U = [? 2((r1+i?h)3 ?h?(??2) + 2(r1+i?h)2 ?(h2)?(??2)]? 8(r1+i?h)4 ?(??2)
? 2(r1+i?h)2 ?(h2)],
V = [ 1h4 + 2((r1+i?h)?h3)]? 8(r1+i?h)4 ?(??2) ? 2(r1+i?h)2 ?(h2)],
W = [ 1((r1+i?h)4 ?(??4)]? 8(r1+i?h)4 ?(??2) ? 2(r1+i?h)2 ?(h2)]
The equation 3.26 must be satisfied at every nodal point of the grid within the
boundary of the plate. When using the Finite Difference Method to solve a partial
differential equation of an annular disc, more accurate results can be obtained by using
higher-order FD approximations than by decreasing the mesh width [8]. Here, a 13 node
finite difference approximation was used in which each node is connected to 13 nodes
other, as shown in figure 3.5. Figure
3.5 Initial and Boundary Conditions
To run a simulation it is first necessary to define the initial and boundary conditions.
It is reasonable to assume that initially the annular plate has no deflection, so the forces
and moments of the plate due to its weight can be neglected [9]. Therefore the initial
conditions of the plate is:
wi,j,k = 0 (3.37)
wi,j,k?1 = 0 (3.38)
26
Table 3.4: Natural frequencies of the disc obtained by the finite difference method
Serial Number Natural Frequency (Hz)
1 979
2 1135
3 1212
4 1827
5 12.98
6 15.831
7 19.1
8 7075
9 7201
It is assumed that the annular disc is clamped at its inner radius, r=r1:
wi,j,k
vextendsinglevextendsingle
vextendsingle(r=r1)
vextendsinglevextendsingle
vextendsingle = 0 (3.39)
?w
?r
vextendsinglevextendsingle
vextendsingle(r=r1)
vextendsinglevextendsingle
vextendsingle = 0 (3.40)
The boundary condition at the outer edge of the disc, i.e at r=r2 is assumed to
be free, so the bending moment and shear force at the outer edge will both be zero.
Mathematically, this can be written as follows:
Vr(r = r2) = ?[D ??r(?2w)+(1??) ?r??(1r ?
2w
?r?? ?
1
r2
?w
?? )] = 0 (3.41)
Mr(r = r2) = ?D[?
2w
?r2 +?(
1
r
?w
?r +
1
r2
?2w
??2 )] = 0 (3.42)
The magnitude of the forcing function is taken to be 5*10(?4) Newton in the static as
well as in the dynamic analysis.
27
Figure 3.6: Response due to a force applied at node 4 of 5X5 grid points
Figure 3.7: Response due to a force applied at node 16 of 5X5 grid points
28
Figure 3.8: Response due to a force applied at node 22 of 5X5 grid points
Figure 3.9: Response due to a force applied at nodes 2,7,12 and 22 of 5X5 grid points
29
Figure 3.10: Response due to a force applied at nodes 3,8,13,18 and 23 of 5X5 grid points
Figure 3.11: Response due to a force applied at nodes 18,19,23 and 24 of 5X5 grid points
30
3.6 Stability of the Algorithm
The stability of the algorithm can be tested by determining whether the iterative
scheme described in equation 3.36 converges to a solution. Discretizing equation 3.26
using equations 3.27-3.35, with zero external force and simplifiying yields:
c[P?wi,j,k +Q?wi+1,j,k +R?wi?1,j,k +S?(wi,j+1,k +wi,j?1,k)+T?(wi+1,j+1,k +wi+1,j?1,k)
+U?(wi?1,j+1,k +wi?1,j?1,k)+V?(wi+2,j,k +wi?2,j,k)+W?(wi,j+2,k +wi,j?2,k)]
= ?wi,j,k+1 +2wi,j,k ?wi,j,k?1
(3.43)
where
P? = [6(??)4 + 6(r1/h+i)4 + 8(??)
2
(r1/h+i)2 ?
8(??)2
(r1/h+i)4 ?
2(??)4
(r1/h+i)2]
(3.44)
Q? = [?4(??)4 ? 4(??)
4
(r1/h+i) +
(??)4
(r1/h+i)3 ?
4(??)2
(r1/h+i)3 ?4
(??)2
(r1/h+i)2]
(3.45)
R? = [?4(??)4 + 4(??)
4
(r1/h+i) ?
(??)4
(r1/h+i)3 +
4(??)2
(r1/h+i)3 ?4
(??)2
(r1/h+i)2]
(3.46)
S? = [? 4(r1/h+i)4 ? 4(??)
2
(r1/h+i)2 +
4(??)2
(r1/h+i)4] (3.47)
31
T? = [ 2(??)
2
(r1/h+i)3 +
2(??)2
(r1/h+i)2] (3.48)
U? = [? 2(??)
2
(r1/h+i)3 +
2(??)2
(r1/h+i)2] (3.49)
V? = [(??)4 + 2(??)
4
(r1/h+i)] (3.50)
W? = [ 1(r1/h+i)4] (3.51)
The deflection wi,j,k can be expressed as:
wi,j,k = ?(t)ei?Rei?? (3.52)
where ? and ? are the spatial frequencies in R and ?, respectively, and i = ??1.
Substituting for wi,j,k from equation 3.52 into equation 3.43 and simplifying yields;
(c?(t))(P?+Q?(ei??R)+R?(e?i??R)+S?(ei??? +e?i???)+T?(ei(??R+??) +ei(??R???))
+U?(ei(???R+??) +e?(i??R+??))+V?(ei2??R +e?i2??R)+W?(ei2??? +e?i2???))
= ??(t+?t)+2?(t)??(t??t)
(3.53)
Using Euler?s identity, ei? = cos(?)+isin(?) in the above equation gives
(c?(t))(P?+Q?(cos(??R))+R?(cos(??R))+S?(2cos(???))+T?(2cos(??R)
cos(???))+U?(2cos(??R)cos(???))+V?(2cos(2??R))+W?(2cos(2???)))
= ??(t+?t)+2?(t)??(t??t)
32
(3.54)
2?c[(P??2V??2W?+Q?(cos(??R))+R?(cos(??R))+S?(2cos(???))+T?(2cos(??R)
cos(???))+U?(2cos(??R)cos(???))+V?(4cos2(??R))+W?(4cos2(???)))]
= ?(t+?t)+?(t??t)?(t)
(3.55)
In general, the solution is assumed to be in the form wi,j,k in equation 3.52 and contains
all multiples of the spatial frequency. Even if these are not present in the initial conditions
or imposed by the boundary conditions, they are likely to be introduced by rounding
errors. This means that the unbounded amplification of ?(t) must be guarded against
all ?[10]. The original error component will not grow in time , provided that
vextendsinglevextendsingle
vextendsingleei(t)
vextendsinglevextendsingle
vextendsingle ? 1
This is a well-known Von neumann stability condition[10, 11]. Therefore, denoting
?(t)=ei(t),the right hand side of the equation can be simplified as:
?(t+?t)+?(t??t)
?(t) =
ei(t+?t) +ei(t??t)
ei(t) = e
i?t +e?i?t
Using Euler?s identity to expand the right hand side of the above equation gives
ei?t = cos(?t)+isin(?t)and
e?i?t = cos(?t)?isin(?t)
33
which leads to
ei?t +e?i?t = 2cos(?t)
The minimum and maximum value of ei?t+e?i?t will be -2 and 2 and cos(?t) will have
minimum and maximum values of -1 and +1, respectively. Since | 2cos(?t) | can only
take a maximum value of 2, the following holds:
?(t+?t)+?(t??t)
?(t) ? 2
Therefore:
| 2?c[(P??2V??2W?)+Q?(cos(??R))+R?(cos(??R))+S?(2cos(???))+T?(2cos(??R)
cos(???))+U?(2cos(??R)cos(???))+V?(4cos2(??R))+W?(4cos2(???))] |? 2
(3.56)
The maximum value each cosine function (cos(??R))(cos(???)),cos2(??R),cos2(???)
can have is 1, when (cos(??R))=(cos(???))= -1. Therefore, the above equation can be
written as
| 2?c[(P??2V??2W?)?Q??R??2S?+2T?+2U?+4V?+4W?] |? 2
(3.57)
34
Figure 3.12: Forcing Function
Substituting for c, P?,Q?,R?,S?,T?,U?,V?andW? into the above equation yields:
| 2?c[16??4 + 16(r
1/h+i)4
+32 ??
2
(r1/h+i)2 ?16
??2
(r1/h+i)4 +2
??4
(r1/h+i)2] |? 2
(3.58)
The sampling time ? t =0.04263 seconds, which satisfies the stability requirement
and is suffcient to cover a broad range of dynamics of the flexible annular plate. The
response of the annular plate is considered over 4.5 seconds.
3.7 Extraction of Eigenvalues and Eigenvectors: Finite Element Approach
Modal analysis can be used to determine the natural frequencies and mode shapes
of a structure. The natural frequencies and modeshapes are important parameters in
the design of a structure for dynamic loading conditions. The eigenvalue and eigenvector
35
Figure 3.13: Dynamic response in transverse direction due to force applied at point 20
of 10X5 grid points
Figure 3.14: Dynamic response in transverse direction due to force applied at point 23
of 10X5 grid points
36
Figure 3.15: Dynamic response in transverse direction due to force applied at point 28
of 10X5 grid points
must therefore be solved for mode-frequency analyses. It takes the following form
[K][?i] = ?i[M][?i] (3.59)
where:
[K]= structure stiffness matrix
[?i]=eigenvector
?i=eigenvalue
[M]= structure mass matrix
ANSYS 8.0 offers several different algorithms to solve the eigenvalue problem like
reduced, block Lanczos, subspace, unsymmetric, damped and QR damped methods.
37
3.7.1 Block Lanczos Method
This method is very commonly used to extract natural frequencies. It uses the
Lanczos algorithm, where the Lanczos recursion is performed for a block of vectors
using a sparse matrix solver. This method is especially powerful when searching for
eigenfrequencies in a defined part of the eigenvalue spectrum of a given system.
3.7.2 Contact Analysis
The contact problem in this study of brake rotor and pad is highly nonlinear. Con-
sequently, it is very important to understand the physics of the problem in order to
run the model more efficiently for accurate results. Two significant problem arose ini-
tially: the region of contact was unknown and the contact problem needed to account
for friction. Friction responses can be chaotic, making solution convergence difficult. In
addition the finite element models of the contacting bodies were generated in such a
way that node-to-node contact was neither achievable nor desirable when contact was
established. For solving this contact type problem, the finite element software ANSYS
makes use of a special type of element CONTAC49 to obtain accurate results. This
element also takes care of the friction between the disc and pad surfaces and can be used
to represent contact and sliding in three dimensions. This element has five nodes and
each node has three degrees of freedom: translations in the nodal x, y and z directions.
This type of element supports elastic coulomb friction and rigid coulomb friction. Both
type of analyses were performed for static as well as dynamic systems. The two po-
tential contact surfaces are referred to as the target surface or the contact surface. The
target surface was represented by target nodes I, J, K and L, and the contact surface
was represented by the contact node M.
38
Table 3.5: Natural frequencies (Hz) of the disc using finite element analysis
Serial Number Natural Frequency,?mn(Hz)
1 998
2 1161
3 1205
4 1840
5 2937
6 4380
7 6212
8 7026
9 7295
Figure 3.16: Another view of the meshed brake rotor model
Figure 3.17: Equivalent disc model
39
Figure 3.18: Modeshape of a disc at natural frequency of 1205 Hz
Figure 3.19: Modeshape of a disc at natural frequency of 1840 Hz
Figure 3.20: Modeshape of a disc at natural frequency of 2937 Hz
40
Figure 3.21: The generation of contact element Contac49 at the disc and pad interface
in ANSYS 8.0
Figure 3.22: Pressurized brake disc model in ANSYS 8.0
41
Chapter 4
Theoretical Modeling of Disc-Brake System
4.1 Lagrange Formulation and Equation of Motion
Although a brake pad and rotor have many modes of vibration, only a single mode
of the complete system contributes to instability and hence noise. In this case only a few
modes of vibration are required to model the system. In this chapter, a simulation model
for the disc-brake system and frictional forces at the disc-pad interface is developed.
Stability analysis has been performed to determine under what parametric condition the
system becomes unstable, assuming that the existence of a limit cycle represents the
noisy state of disc brake system [24]. In order to study the nonlinear behavior of the
brake system and its stability, a mathematical model with the pad and disk modeled as
an idealized mass-spring-damper system was constructed. Considering the disc-rotor to
be an annular disc, it can be characterized by a dynamic system of differential equations
of order four. Therefore, a mathematical model with 4 degrees of freedom is presented.
The degrees of freedom x1, y1, x2, y2 are the inplane and transverse displacements of the
pad and disk, respectively, in the cartesian co-ordinate system as shown in figure 4.1.
This is a minimal DOF model that includes both in-plane and transverse vibrations of the
brake components in order to obtain phase portraits and time histories in a reasonable
amount of time. The resilience and the dissipation of the disk and the pad materials are
modeled by linear springs and dampers. The linear springs and dampers for the pad are
denoted by k1 and C1, respectively, and same for the disc by k2 and C2, respectively.
42
The contact between the disc and pad is modeled by a transverse spring and damper
denoted by kT and CT as shown in figure 4.1.
In order to derive the equations of motion, the kinetic energy T, the potential energy
V and the dissipative energy D of the system, can be written as:
T = 12m1x21 + 12m1y21 + 12m2x22 + 12m2y22 (4.1)
V = 12k1x21 + 12k1x22 + 12kT(y1 ?y2)2 (4.2)
D = 12C1x?21 + 12C2x?22 + 12CT(y?1 ?y?2)2 (4.3)
where m1 and m2 are the masses of the pad and disc respectively.
The generalized co-ordinates are taken to be q1= x1, q2=y1 for the pad and q3=x2,
q4=y2 for the disc.
The general Lagrange formulation used to derive equations of motion with respect
to generalized co-ordinates is:
d
dt(
?T
?q?i)+
?D
?q?i +
?V
?qi = Qi, (4.4)
where Qi is the generalized force associated to the ith generalized coordinate. The
generalized forces can thus be expressed as : Q1 = Ff, Q2 = Nstatic, Q3 = -Ff, Q4 =0,
where Ff is the friction force and Nstatic is the initial force acting on the system.
Finally, the equations of motion can be written as:
m1x?1 +C1x?1 +k1x1 = Ff (4.5)
43
Figure 4.1: Four degrees of freedom model
m1y?1 +CT(y?1 ?y?2)+kT(y1 ?y2)+kT(y1 ?y2)3 = Nstatic (4.6)
m2x?2 +C2x?2 +k2x2 = ?Ff (4.7)
m2y?2 ?CT(y?1 ?y?2)?kT(y1 ?y2)?kT(y1 ?y2)3 = 0 (4.8)
The cubic non-linearity is added in the above two equation with reference to Jear-
siripongkul and Hagedorn.[23].
During braking the system experiences both in-plane and out-of-plane vibrations. To
take into account the transverse vibration of the system, the normal force acting on the
system is considered to be [12] :
N = Nstatic +Ndynamic (4.9)
44
where the static normal force depends on the pressure applied by the piston on the pad,
Nstatic = P ? S and the dynamic component is modeled as:
Ndynamic = kT(y1 ?y2)+CT(y?1 ?y?2) (4.10)
4.1.1 Friction Models
To analyze the behavior of the two systems two friction models were considered.
The first was a Coulomb type friction force independent of the magnitude of velocity
and always opposing the relative motion. Suppose the friction force comes from the
dynamic coefficient (the relative movement between pad and disc, expressed by the
relative velocity vr negationslash= 0), then one can write:
Ff = ??d(Nstatic +Ndynamic)sgn(vr) (4.11)
A more realistic model for the friction takes into consideration the variation of the
dynamic coefficient with the relative velocity using the formula ?d=?s - ?vr. In the
stick phase, as the relative velocity (vr=0) between the pad and disc approaches zero,
the friction force is modeled as the minimum value between the in-plane applied force FT
and the static friction force Fs=?s Nstatic. It cannot be bigger than the applied in-plane
force FT and must oppose it.
Considering the relative motion between pad and disc, the applied in-plane force is
evaluated as [24]
45
FT = ?k1x1 ?c1x?1 +k2x2 +c2x?2 (4.12)
In the slip phase, where the relative movement is significant (vr negationslash= 0), the friction force
is ?dN sgn(vr), a coulomb type friction force. Finally, the stick-slip friction model that
induces excitation in the system is
if vr=0, stick condition
Ff=-min(| FT,Fs |)sgn(FT)
and if, vr negationslash= 0, slip condition
Ff=(?s - ? vr)(Nstatic+Ndynamic)sgn(vr)
4.2 Numerical Simulation and Stability Analysis
The equations of motion are coupled differential equations of second order:
r?= F(t,r,s,r?,s?)
s?= G(t,r,s,r?,s?)
where?denotes differentiation with respect to time.
Assuming
r?= u
s?= v
46
Therefore
u?= F(t,r,s,u,v)
v?= G(t,r,s,u,v)
reducing the second order differential equations to a system of four simultaneous equa-
tions of first order. A single integration step, h of time t, is then given by,
rn+1 = rn +(1/6)?k1 +(1/3)?k2 +(1/3)?k3 +(1/6)?k4 +O(h5)
sn+1 = sn +(1/6)?l1 +(1/3)?l2 +(1/3)?l3 +(1/6)?l4 +O(h5)
un+1 = un +(1/6)?m1 +(1/3)?m2 +(1/3)?m3 +(1/6)?m4 +O(h5)
vn+1 = vn +(1/6)?p1 +(1/3)?p2 +(1/3)?p3 +(1/6)?p4 +O(h5)
where
k1 = h?u (4.13)
k2 = h?(u+m1/2)
k3 = h?(u+m2/2)
k4 = h?(u+m3)
l1 = h?v
l2 = h?(v +p1/2)
l3 = h?(v +p2/2)
l4 = h?(v +p3)
47
Figure 4.2: (a)Motion & (b)Time history of the disc (Force N=1 and b=0.1)respectively
m1 = h?F(t,r,s,u,v)
m2 = h?F(t+h/2,r +k1/2,s+l1/2,u+m1/2,v +p1/2)
m3 = h?F(t+h/2,r +k2/2,s+l2/2,u+m2/2,v +p2/2)
m4 = h?F(t+h,r +k3,s+l3/2,u+m3,v +p3/2)
p1 = h?G(t,r,s,u,v)
p2 = h?G(t+h/2,r +k1/2,s+l1/2,u+m1/2,v +p1/2)
p3 = h?G(t+h/2,r +k2/2,s+l2/2,u+m2/2,v +p2/2)
p4 = h?G(t+h,r +k3,s+l3,u+m3,v +p3)
The above method to solve the coupled differential equation is known as Runge-Kutta?s
JBM method. A C-program can be written to run the simulation and various parameters
are used to predict the behavior. For this study, various simulations were run to obtain
the stability and instability conditions as shown in figures 4.2 to 4.12. A parametric
study was done using different values of static coefficient of friction, damping , initial
velocity, dynamic coefficient of friction, force magnitudes and damping of disc & pad. In
figures 4.2-4.3 the various parameters used are mass m1=m2=1, damping c1=c2=0.1 &
cT=0.5, stiffness k1=k2=kT=1, initial velocity u=1, force magnitude N=1 and ?s=0.01
and for various values of ?. In figures 4.4-4.5 the various parameters used are mass
48
Figure 4.3: (a)Motion & (b)Time history of the disc (Force N=1 and b=1)respectively
Figure 4.4: (a)Motion & (b)Time history of the disc (Force N=1 and
us=0.01)respectively
m1=m2=1 damping c1=c2=0.1 & cT=0.5, stiffness k1=k2=kT=1, initial velocity u=1,
force magnitude N=1, ?=0.01 and various values of ?s.
Figure 4.5: (a)Motion & (b)Time history of the disc respectively (Force N=1 and us=0.1)
49
Chapter 5
Results and Conclusions
5.1 Summary of Results
5.1.1 Analytical Method
The natural frequencies obtained through the analytical method were in generally
good with the natural frequencies obtained through other methods given in table 5.1.
Table 5.2 shows the percentage differences of the natural frequencies of the various
Table 5.1: Natural frequencies of the disc obtained through all methods
Serial Number Analytical Finite Difference Finite Element Experimental
1 979 952 998 963
2 1135 1152 1161 1184
3 1212 1198 1205 1235
4 1827 1879 1840 1901
5 2902 3012 2937 2933
6 4316 4287 4380 4358
7 6282 6159 6212 6104
8 7075 6995 7026 7098
9 7201 7197 7295 7310
methods with respect to the analytical method.
5.2 Finite Difference Method
A numerical method of solution of the governing partial differential equation de-
scribing the characteristic behavior of a flexible annular plate has been presented and
discussed. A simulation algorithm characterizing the behavior of a flexible annular plate
50
Table 5.2: Percentage deviation of natural frequencies of the disc with respect to the
experimental method
Serial Number Analytical Finite Difference Finite Element
1 -1.661 1.142 -3.634
2 4.138 2.702 1.942
3 1.862 2.995 2.429
4 3.892 1.157 3.208
5 1.056 -2.693 -0.136
6 0.963 1.629 -0.504
7 -2.916 -0.901 -1.769
8 0.324 1.451 1.014
9 1.491 1.545 0.205
with inner edge clamped or fixed and outer edges free has been developed in C pro-
gramming and Matlab environment. The simulation results has been presented in time
domain. The deflection of annular disc has been studied by varying the number of sec-
tions, n in the radial direction and m in the circumferential direction and the results
have been validated by comparing these with previously reported results using other
techniques. The stability of the algorithm has been analyzed and sufficient condition for
stability of the algorithm has been obtained.
5.2.1 Case I
Figure 5.1 indicates the response in the form of displacement. Various mesh grids
generated using the finite difference method and the disc model were simulated for dif-
ferent constant forcing functions. An example of 5X5 grid points and a constant force of
magnitude 5*10(?4)Newton applied at node 17 of the disc is shown here. It was found
that the maximum reponse occurs at this point on the disc. The main advantage of this
method is that the response can be found at any point of the disc.
51
Figure 5.1: Response due to a force applied at node 17 of 5X5 grid points
5.2.2 Case II
The plot shown in figure 5.3 indicates the response in the form of displacement.
Various mesh grids were generated using the central difference method and the disc
model was simulated for different time varying forces. The example of 10X5 grid points
and a time varying forces of magnitude 5*10(?4) newton applied at node 17 of the disc
is shown here. The plot of the forcing function is shown in figure 5.2 and the response is
shown below in figure 5.3. The advantage of this method is that any practically possible
time varying forcing function can be applied to the disc and the response determined.
In the above two cases using the finite difference method, friction was not taken into
account.
5.3 Finite Element Method
Friction was, however, taken into account for the contact analysis. A specific type of
element CONTAC49, was used to conduct this analysis. This element supports friction
52
Figure 5.2: Forcing Function
Figure 5.3: Dynamic response in transverse direction due to a force applied at point 17
of 10X5 grid points
53
Figure 5.4: Nodal displacement of the disc in the z direction
and the results obtained are shown in the figure 5.4 below. The light blue color in the
figure 5.4 shows the maximum displacement of disc in z direction and figure 5.5 shows
the reponse of disc in the z direction with respect to time. It was found that all of the
disc nodes which were in contact with the pad had almost the same displacement in the
z direction.
5.4 Four Degrees of Freedom Model
A four degree-of-freedom model was used to investigate the basic mechanisms of an
instability that is one of the causes of disc brake noise. The damping of disc and pad
are of equally importance as it reduces the chances of instability. The model was also
used to determine the conditions necessary to prevent this instability. The numerical
analysis suggests that when the natural frequency of disc and pad are close then a brake
is more likely to be noisy. A non-linear analysis of the coupled equation was conducted to
investigate the dynamical behavior of model for various system and friction parameters.
54
Figure 5.5: Reponse in the z direction with repect to time
The numerical analysis suggested that the initial velocity of the disc also plays a major
role in the instability of the disc brake system.
5.5 Conclusions
An analytical method was developed and is reported here that can be used to predict
the natural frequencies of an annular disc with boundary conditions. A new methodol-
ogy has been presented for studying the non-linear dynamics of the annular disc. An
important advantage of the new methodology is that it does not require the often dif-
ficult task of inducing squeal in the laboratory. The analysis presented using the finite
difference method is independent of the laws of friction.
Finite element analysis was used to incorporate the sliding friction between the rotor
and pad surfaces and the dynamic response obtained through this analysis has impor-
tant implications. The results obtained through this analysis identify the locations of
maximum deflection in a pressurized brake system. The reduction of these deflections
55
either partially or completely can thus eliminate the instability in a disc brake system.
An analytical four degrees of freedom model has been used to identify the necessary
conditions for stability. The initial velocity, the dynamic coefficient of friction, damping
of the disc and pad and the magnitude of the force all play major roles in the stability
of the disc brake system.
56
Bibliography
[1] Hagedorn,P., ?Modeling Disk Brakes with respect to squeal?, TU Darmstadt, Ger-
many
[2] Akay A., ?Acoustics of friction?, Journal of the Acoustic Society of America 2002,
111(4):1525-1548.
[3] Dunlap, K.B.,Riehle, M.A. and Longhouse, R.E., ?An investigative overview of
automotive disk brake noise?, SAE Paper 1999-01-0142.
[4] Cao,Q., Ouyang,H., Friswell, M.I. and Mottershead, J.E., ?Linear eigenvalue analy-
sis of the disc-brake squeal?, International Journal For Numerical Methods In En-
gineering 2004.
[5] Fischer, M. and Bendel, K., ?Extensive 3-D Vibrometry on Brake Disks?,Corporate
Research and Development Applied Physics,Polytec LM Info special issue 2004.
[6] Mote, C.D., ?Stability of Circular Plates Subjected to Moving Loads?, J.Franklin
Inst., 290, pp.329-344, 1970.
[7] North, N.R., ? Disc Brake Squeal? Proceedings of IMechE 1976;(C38/76):169-176.
[8] Abramowithz, M., and Cahill, W.F., ?On the vibration of the square clamped
plate?,Journal of Association for Computing Machinery, 2(3), pp.162-168, 1995.
[9] Timoshenko, S. and Woinowsky-Krieger, S. ,?Theory of Plate and Shells?, McGraw-
Hill Book Company, New York.
[10] Azad, A.K.M., ?Analysis and design of control mechanisms for flexible manipu-
lator systems?,PhD.Thesis Dept. of Automatic Control and Systems Engineering,
University of Sheffield, UK, 1994.
[11] Vemuri,M.O. and Karplus, W.J., ?Digital computer treatment of partial differential
equations?,Prentice-Hall Series in Computational Mathematics, New Jersey, 1994.
[12] Chung,C.H. and Kobayashi, K., ?A new analysis method for brake squeal part 1:
Theory for Modal Domain Formulation and Stability Analysis?, Society of Auto-
motive Engineers,Inc., 2000.
[13] Yuan, Y.A., ?A study of the effects of negative friction-speed slope on brake squeal?,
In Proceedings of 1995 ASME Design Engineering Conference, Boston,1995, Vol.3,
Part A, 1135-1162.
57
[14] Mills H.R., ?Brake Squeak?, Institute of Automobile Engineers, Report 9000B, 1938.
[15] Ouyang,H., Mottershead, J.E., Cartmell, M.P. and Brookfield, D.J., ?Friction-
induced vibration of an elastic slider on a vibrating disc?, Int. J. Mech. Sct., 1999,
41, 325-336.
[16] Ouyang,H., Cao Q., Mottershead, J.E. and Treyde, T., ?Vibration and squeal of a
disc brake: Modeling and experimental results?, Proc. Instn. Mech. Engrs Vol 217
Part D:J.Automobile Engineering, 867-875.
[17] Ouyang,H. and Mottershead, J.E., ?Disk under the action of a rotating friction
couple?, Journal of Applied Mechanics Vol 71 November 2004, 753-758.
[18] Chowdhary, Harsh V., Bajaj, K. Anil, Krousgrill, Charles M., ?An Analytical
Approach To Model Disc Brake System for Squeal Prediction?, Proceedings of
DETC?01.
[19] McDaniel, J., James, M., Chen and Clarke, C. L., ?Acoustic radiation Models of
Brake Systems from Stationary LDV Measurements?, Proceedings of IMEC 99,
USA.
[20] Reismann H., ?Elastic Plates: Theory and Application?, John Wiley and Sons,
1988.
[21] Young, T.H., Lin,C.Y., ?Stability of a spinning disk under a stationary oscillat-
ing unit?, Department of Mechanical Engineering, National Taiwan University of
Science and Technology, Taipei, 106, Taiwan.
[22] Dav, G., and Milazzo, A., ?A new symmetric and positive definite boundary element
formulation for lateral vibration of plates?, Journal of Sound and Vibrations,Vol.
206, No.4, 507-521, 1997.
[23] Jearsiripongkul,T., Hagedorn, P.,?Parametric Stability Analysis of a floating Caliper
Disk Brake?,Institut fu?r Mechanik, TU-Darmstadt.
[24] Shin, K., Brennan, M.J., Oh, J.E. and Harris,C.J., 2002 ,?Analysis of disk brake
noise using a two-degrees-of-freedom model?, Journal of Sound and Vibration,
Vol.254(5), pp.837-848.
58
Appendices
59
Appendix A
C-program to solve nonlinear differential equation using finite
difference method
Program to solve the differential equation
include< stdio.h >
include< math.h >
include< stdlib.h >
define pi 3.147
define v 10
define w 5
int jadeja(int M, int N, int p, int q,double x,double r1,double r2,int **stcl,
double **val);
int main()
( FILE ?vijay1,?vijaya;
int inx,p=0,q=0,z=0,r=0,b=0;
int ix,jx,ix,jx,i,j,xp,y,yp;
int k=0;
int l=0;
int M=0; //Number of divisions//
int N=0; //Number of circles//
double r1,r2,x;
int mnx=0,nnx=0,jy = 0;
int inx = 0,jnx = 0;
int px=0,qx=0, row = 0, col = 0;
vijay1=fopen(?h://maharaja//raheja.txt?,?r?);
vijaya=fopen(?h://maharaja//vishu.txt?,?w?);
fscanf(vijay1,?d/n/d/n/lf/nlf/n?,&M,&N,&r1,&r2);
printf(?M=d/tN=d/tr1=lf/tr2=lf/n?,M,N,r1,r2);
int *stcl=NULL;
double *val=NULL; double **anx=NULL;
x=(r2-r1)/M;
anx=(double**)malloc((129)*sizeof(double*));
for(inx = 1;inx <= 129;inx ++)
(
anx[inx]=(double*)malloc((129)*sizeof(double));
)
for(r=1;r<=50;r++)
60
(
for(y=1;y<=50;y++)
(
anx[r][y]=0.0;
)
)
for(p=1;p<=N;p++)
(
for(q=1;q<=M;q++)
(
free (stcl);
free (val);
val = NULL;
stcl = NULL;
jadeja(M,N,p,q,x,r1,r2,&stcl,&val);
row=(p-1)*M+q;
for ( i=1;i<27;i=i+2 )
(
col=(stcl[i]-1)*M+stcl[i+1];
anx[row][col] =+ val[(i/2)+1];
if(q==1)
(
printf(?lf/tnode[d][d]/t?,val[(i/2)+1],stcl[i],stcl[i+1]);
anx[row][col]=0;
)
else if(q==2 && anx[row][col]==val[1])
(
anx[row][col]=val[1]+val[7];
)
else if(q==M-1 && anx[row][col]==val[1])
(
anx[row][col]=val[1]+val[18];
)
else if(q==M-1 && anx[row][col]==val[8])
(
anx[row][col]=val[8]+val[15];
)
else if(q==M-1 && anx[row][col]==val[10])
(
anx[row][col]=val[10]+val[16];
)
else if(q==M-1 && anx[row][col]==val[13])
61
(
anx[row][col]=val[13]+val[17];
)
else if(q==M && anx[row][col]==val[1])
(
anx[row][col]=+val[15];
)
else if(q==M && anx[row][col]==val[2])
(
anx[row][col]=+val[16];
)
else if(q==M && anx[row][col]==val[4])
(
anx[row][col]=+val[17];
)
else if(q==M && anx[row][col]==val[6])
(
anx[row][col]=+val[18];
)
)
for(r=1;r<=M*N;r++)
(
for(y=1;y<=M*N;y++)
(
printf(?lf/t?,anx[r][y]);
)
printf(?/n?);
)
)
)
fclose(vijay1);
fclose(vijaya);
)
int jadeja(int M, int N, int p,int q,double x,double r1,double r2,
int **stcl,double **val)
(
int i,j,k,z;
int ip;
val=(double*)malloc((29)*sizeof(double));
int I=0;
double dta=(2*pi/N);
62
double rho=7800,deltaT=0.1,c=0,omega=1,mu=0.3;
stcl=(int*)malloc((26)*sizeof(int));
if(p==1 && q==1)
(
(*stcl)[1]=p;(*stcl)[2]=q;(*stcl)[3]=p+1;(*stcl)[4]=q;(*stcl)[5]=p+2;
(*stcl)[6]=q;(*stcl)[7]=N;(*stcl)[8]=q;(*stcl)[9]=N-1;(*stcl)[10]=q;
(*stcl)[11]=-10;(*stcl)[12]=-10;(*stcl)[13]=-10;(*stcl)[14]=-10;
(*stcl)[15]=p;(*stcl)[16]=q+1;(*stcl)[17]=p;(*stcl)[18]=q+2;
(*stcl)[19]=p+1;(*stcl)[20]=q+1;(*stcl)[21]=-10;(*stcl)[22]=-10;
(*stcl)[23]=-10; (*stcl)[24]=-10; (*stcl)[25]=N; (*stcl)[26]=q+1;
)
else if(p==1 && q==2)
(
(*stcl)[1]=p;(*stcl)[2]=q;(*stcl)[3]=p+1;(*stcl)[4]=q;(*stcl)[5]=p+2;
(*stcl)[6]=q;(*stcl)[7]=N;(*stcl)[8]=q;(*stcl)[9]=N-1;(*stcl)[10]=q;
(*stcl)[11]=p;(*stcl)[12]=q-1;(*stcl)[13]=-1;(*stcl)[14]=-1;
(*stcl)[15]=p;(*stcl)[16]=q+1;(*stcl)[17]=p;(*stcl)[18]=q+2;
(*stcl)[19]=p+1;(*stcl)[20]=q+1;(*stcl)[21]=p+1;(*stcl)[22]=q-1;
(*stcl)[23]=N; (*stcl)[24]=q-1; (*stcl)[25]=N; (*stcl)[26]=q+1;
)
else if( p==N && q==1)
(
(*stcl)[1]=p;(*stcl)[2]=q;(*stcl)[3]=1;(*stcl)[4]=q;(*stcl)[5]=2;
(*stcl)[6]=q;(*stcl)[7]=p-1;(*stcl)[8]=q;(*stcl)[9]=p-2;(*stcl)[10]=q;
(*stcl)[11]=-1;(*stcl)[12]=-1;(*stcl)[13]=-1;(*stcl)[14]=-1;(*stcl)[15]=p;
(*stcl)[16]=q+1;(*stcl)[17]=p;(*stcl)[18]=q+2;(*stcl)[19]=1;(*stcl)[20]=q+1;
(*stcl)[21]=-1;(*stcl)[22]=-1;(*stcl)[23]=-1;(*stcl)[24]=-1;(*stcl)[25]=p-1;
(*stcl)[26]=q+1;
)
else if( p==N-1 && q==1)
(
(*stcl)[1]=p;(*stcl)[2]=q;(*stcl)[3]=p+1;(*stcl)[4]=q;(*stcl)[5]=1;
(*stcl)[6]=q;(*stcl)[7]=p-1;(*stcl)[8]=q;(*stcl)[9]=p-2;(*stcl)[10]=q;
(*stcl)[11]=-1;(*stcl)[12]=-1;(*stcl)[13]=-1;(*stcl)[14]=-1;(*stcl)[15]=p;
(*stcl)[16]=q+1;(*stcl)[17]=p;(*stcl)[18]=q+2;(*stcl)[19]=N;(*stcl)[20]=q+1;
(*stcl)[21]=-1;(*stcl)[22]=-1;(*stcl)[23]=-1;(*stcl)[24]=-1;(*stcl)[25]=p-1;
(*stcl)[26]=q+1;
)
else if(p==1 && q>2)
(
if(q>2 && q ?=(M-2))
(
63
(*stcl)[1]=p;(*stcl)[2]=q;(*stcl)[3]=p+1;(*stcl)[4]=q;(*stcl)[5]=p+2;
(*stcl)[6]=q;(*stcl)[7]=N;(*stcl)[8]=q; (*stcl)[9]=N-1;(*stcl)[10]=q;
(*stcl)[11]=p;(*stcl)[12]=q-1;(*stcl)[13]=p;(*stcl)[14]=q-2;(*stcl)[15]=p;
(*stcl)[16]=q+1;(*stcl)[17]=p; (*stcl)[18]=q+2;(*stcl)[19]=p+1;
(*stcl)[20]=q+1;(*stcl)[21]=p+1;(*stcl)[22]=q-1;(*stcl)[23]=N;
(*stcl)[24]=q-1;(*stcl)[25]=N; (*stcl)[26]=q+1;
)
else if(p==1 && q==M)
(
(*stcl)[1]=p;(*stcl)[2]=q;(*stcl)[3]=p+1;(*stcl)[4]=q;(*stcl)[5]=p+2;
(*stcl)[6]=q;(*stcl)[7]=N;(*stcl)[8]=q; (*stcl)[9]=N-1;(*stcl)[10]=q;
(*stcl)[11]=p;(*stcl)[12]=q-1;(*stcl)[13]=p;(*stcl)[14]=q-2;(*stcl)[15]=-1;
(*stcl)[16]=-1;(*stcl)[17]=-2;(*stcl)[18]=-2;(*stcl)[19]=-1;(*stcl)[20]=-1;
(*stcl)[21]=p+1;(*stcl)[22]=q-1;(*stcl)[23]=N;(*stcl)[24]=q-1;(*stcl)[25]=-1;
(*stcl)[26]=-1;
)
else if (p==1 && q==(M-1))
(
(*stcl)[1]=p;(*stcl)[2]=q;(*stcl)[3]=p+1;(*stcl)[4]=q;(*stcl)[5]=p+2;
(*stcl)[6]=q;(*stcl)[7]=N;(*stcl)[8]=q; (*stcl)[9]=N-1;(*stcl)[10]=q;
(*stcl)[11]=p;(*stcl)[12]=q-1;(*stcl)[13]=p;(*stcl)[14]=q-2;(*stcl)[15]=p;
(*stcl)[16]=q+1;(*stcl)[17]=-1;(*stcl)[18]=-1;(*stcl)[19]=p+1;(*stcl)[20]=q+1;
(*stcl)[21]=p+1;(*stcl)[22]=q-1;(*stcl)[23]=N;(*stcl)[24]=q-1;(*stcl)[25]=N;
(*stcl)[26]=q+1;
)
)
else if(p>=1 && q==M)
(
if(p==2 && q==M)
(
(*stcl)[1]=p;(*stcl)[2]=q;(*stcl)[3]=p+1;(*stcl)[4]=q;(*stcl)[5]=p+2;(*stcl)[6]=q;
(*stcl)[7]=p-1;(*stcl)[8]=q; (*stcl)[9]=N;(*stcl)[10]=q;(*stcl)[11]=p;(*stcl)[12]=q-1;
(*stcl)[13]=p;(*stcl)[14]=q-2;(*stcl)[15]=-1;(*stcl)[16]=-1; (*stcl)[17]=-1; (*stcl)[18]=-1;
(*stcl)[19]=-2;(*stcl)[20]=-2;(*stcl)[21]=p+1;(*stcl)[22]=q-1;(*stcl)[23]=p-1;
(*stcl)[24]=q-1;(*stcl)[25]=-1; (*stcl)[26]=-1;
)
else if(p>2 && p<=(N-2))
(
(*stcl)[1]=p;(*stcl)[2]=q;(*stcl)[3]=p+1;(*stcl)[4]=q;(*stcl)[5]=p+2;(*stcl)[6]=q;
(*stcl)[7]=p-1;(*stcl)[8]=q; (*stcl)[9]=p-2;(*stcl)[10]=q;(*stcl)[11]=p;(*stcl)[12]=q-1;
(*stcl)[13]=p;(*stcl)[14]=q-2;(*stcl)[15]=-1;(*stcl)[16]=-1;(*stcl)[17]=-1;(*stcl)[18]=-1;
(*stcl)[19]=-2;(*stcl)[20]=-2; (*stcl)[21]=p+1;(*stcl)[22]=q-1;(*stcl)[23]=p-1;
64
(*stcl)[24]=q-1;(*stcl)[25]=-1;(*stcl)[26]=-1;
)
else if(p==(N-1) && q==M)
(
(*stcl)[1]=p;(*stcl)[2]=q;(*stcl)[3]=p+1;(*stcl)[4]=q;(*stcl)[5]=1;(*stcl)[6]=q;
(*stcl)[7]=p-1;(*stcl)[8]=q; (*stcl)[9]=p-2;(*stcl)[10]=q;(*stcl)[11]=p;(*stcl)[12]=q-1;
(*stcl)[13]=p;(*stcl)[14]=q-2;(*stcl)[15]=-1;(*stcl)[16]=-1;(*stcl)[17]=-1;(*stcl)[18]=-1;
(*stcl)[19]=-1;(*stcl)[20]=-1;(*stcl)[21]=p+1;(*stcl)[22]=q-1;(*stcl)[23]=p-1;
(*stcl)[24]=q-1;(*stcl)[25]=-1; (*stcl)[26]=-1;
)
else if(p==N && q==M)
(
(*stcl)[1]=p;(*stcl)[2]=q;(*stcl)[3]=1;(*stcl)[4]=q;(*stcl)[5]=2;(*stcl)[6]=q;
(*stcl)[7]=p-1;(*stcl)[8]=q; (*stcl)[9]=p-2;(*stcl)[10]=q;(*stcl)[11]=p;(*stcl)[12]=q-1;
(*stcl)[13]=p;(*stcl)[14]=q-2;(*stcl)[15]=-1;(*stcl)[16]=-1;(*stcl)[17]=-1;
(*stcl)[18]=-1;(*stcl)[19]=-1;(*stcl)[20]=-1;(*stcl)[21]=1;(*stcl)[22]=q-1;
(*stcl)[23]=p-1;(*stcl)[24]=q-1;(*stcl)[25]=-1; (*stcl)[26]=-1;
)
)
else if(p>=2 && q>=1 && q<(M-2))
(
if(p?2 && p?(N-2) && q?2)
(
(*stcl)[1]=p;(*stcl)[2]=q;(*stcl)[3]=p+1;(*stcl)[4]=q;(*stcl)[5]=p+2;(*stcl)[6]=q;
(*stcl)[7]=p-1;(*stcl)[8]=q; (*stcl)[9]=p-2;(*stcl)[10]=q;(*stcl)[11]=p;(*stcl)[12]=q-1;
(*stcl)[13]=p;(*stcl)[14]=q-2;(*stcl)[15]=p;(*stcl)[16]=q+1;(*stcl)[17]=p;(*stcl)[18]=q+2;
(*stcl)[19]=p+1;(*stcl)[20]=q+1;(*stcl)[21]=p+1;(*stcl)[22]=q-1;(*stcl)[23]=p-1;
(*stcl)[24]=q-1; (*stcl)[25]=p-1; (*stcl)[26]=q+1;
)
else if(p==2 && q==2)
(
(*stcl)[1]=p;(*stcl)[2]=q;(*stcl)[3]=p+1;(*stcl)[4]=q;(*stcl)[5]=p+2;(*stcl)[6]=q;
(*stcl)[7]=p-1;(*stcl)[8]=q; (*stcl)[9]=N;(*stcl)[10]=q;(*stcl)[11]=p;(*stcl)[12]=q-1;
(*stcl)[13]=-1;(*stcl)[14]=-1;(*stcl)[15]=p;(*stcl)[16]=q+1;(*stcl)[17]=p;(*stcl)[18]=q+2;
(*stcl)[19]=p+1;(*stcl)[20]=q+1;(*stcl)[21]=p+1;(*stcl)[22]=q-1;(*stcl)[23]=p-1;
(*stcl)[24]=q-1; (*stcl)[25]=p-1; (*stcl)[26]=q+1;
)
else if(p>2 && p<=(N-2) && q==2)
(
(*stcl)[1]=p;(*stcl)[2]=q;(*stcl)[3]=p+1;(*stcl)[4]=q;(*stcl)[5]=p+2;(*stcl)[6]=q;
(*stcl)[7]=p-1;(*stcl)[8]=q; (*stcl)[9]=p-2;(*stcl)[10]=q;(*stcl)[11]=p;(*stcl)[12]=q-1;
(*stcl)[13]=-1;(*stcl)[14]=-1;(*stcl)[15]=p;(*stcl)[16]=q+1;(*stcl)[17]=p;(*stcl)[18]=q+2;
65
(*stcl)[19]=p+1;(*stcl)[20]=q+1;(*stcl)[21]=p+1;(*stcl)[22]=q-1;(*stcl)[23]=p-1;
(*stcl)[24]=q-1;(*stcl)[25]=p-1; (*stcl)[26]=q+1;
)
else if(p==N && q==2)
(
(*stcl)[1]=p;(*stcl)[2]=q;(*stcl)[3]=1;(*stcl)[4]=q;(*stcl)[5]=2;(*stcl)[6]=q;
(*stcl)[7]=p-1;(*stcl)[8]=q; (*stcl)[9]=p-2;(*stcl)[10]=q;(*stcl)[11]=p;(*stcl)[12]=q-1;
(*stcl)[13]=-1;(*stcl)[14]=-1;(*stcl)[15]=p;(*stcl)[16]=q+1;(*stcl)[17]=p;(*stcl)[18]=q+2;
(*stcl)[19]=1;(*stcl)[20]=q+1;(*stcl)[21]=1;(*stcl)[22]=q-1;(*stcl)[23]=p-1;
(*stcl)[24]=q-1;(*stcl)[25]=p-1; (*stcl)[26]=q+1;
)
else if(p==2 && q==1)
(
(*stcl)[1]=p;(*stcl)[2]=q;(*stcl)[3]=p+1;(*stcl)[4]=q;(*stcl)[5]=p+2;(*stcl)[6]=q;
(*stcl)[7]=p-1;(*stcl)[8]=q; (*stcl)[9]=N;(*stcl)[10]=q;(*stcl)[11]=-1;(*stcl)[12]=-1;
(*stcl)[13]=-1;(*stcl)[14]=-1;(*stcl)[15]=p;(*stcl)[16]=q+1;(*stcl)[17]=p;(*stcl)[18]=q+2;
(*stcl)[19]=p+1;(*stcl)[20]=q+1;(*stcl)[21]=-1;(*stcl)[22]=-1;(*stcl)[23]=-1;
(*stcl)[24]=-1;(*stcl)[25]=p-1; (*stcl)[26]=q+1;
)
else if(p>2 && p<=(N-2) && q==1)
(
(*stcl)[1]=p;(*stcl)[2]=q;(*stcl)[3]=p+1;(*stcl)[4]=q;(*stcl)[5]=p+2;(*stcl)[6]=q;
(*stcl)[7]=p-1;(*stcl)[8]=q; (*stcl)[9]=p-2;(*stcl)[10]=q;(*stcl)[11]=-1;(*stcl)[12]=-1;
(*stcl)[13]=-1;(*stcl)[14]=-1;(*stcl)[15]=p;(*stcl)[16]=q+1;(*stcl)[17]=p;
(*stcl)[18]=q+2;(*stcl)[19]=p+1;(*stcl)[20]=q+1;(*stcl)[21]=-1;(*stcl)[22]=-1;
(*stcl)[23]=-1;(*stcl)[24]=-1;(*stcl)[25]=p-1;(*stcl)[26]=q+1;
)
else if (p<=(N-2) && q>2 && q<M-2)
(
(*stcl)[1]=p;(*stcl)[2]=q;(*stcl)[3]=p+1;(*stcl)[4]=q;(*stcl)[5]=p+2;(*stcl)[6]=q;
(*stcl)[7]=p-1;(*stcl)[8]=q; (*stcl)[9]=p-2;(*stcl)[10]=q;(*stcl)[11]=p;(*stcl)[12]=q-1;
(*stcl)[13]=p;(*stcl)[14]=q-2;(*stcl)[15]=p;(*stcl)[16]=q+1;(*stcl)[17]=p;(*stcl)[18]=q+2;
(*stcl)[19]=p+1;(*stcl)[20]=q+1;(*stcl)[21]=p+1;(*stcl)[22]=q-1;(*stcl)[23]=p-1;
(*stcl)[24]=q-1;(*stcl)[25]=p-1; (*stcl)[26]=q+1;
)
)
else if(p>=2 && p<=(N-2) && q>2 && q<=(M-2))
(
if(p==2)
(
(*stcl)[1]=p;(*stcl)[2]=q;(*stcl)[3]=p+1;(*stcl)[4]=q;(*stcl)[5]=p+2;(*stcl)[6]=q;
(*stcl)[7]=p-1;(*stcl)[8]=q; (*stcl)[9]=N;(*stcl)[10]=q;(*stcl)[11]=p;(*stcl)[12]=q-1;
66
(*stcl)[13]=p;(*stcl)[14]=q-2;(*stcl)[15]=p;(*stcl)[16]=q+1;(*stcl)[17]=p;(*stcl)[18]=q+2;
(*stcl)[19]=p+1;(*stcl)[20]=q+1;(*stcl)[21]=p+1;(*stcl)[22]=q-1;(*stcl)[23]=p-1;
(*stcl)[24]=q-1;(*stcl)[25]=p-1;(*stcl)[26]=q+1;
)
else if(p>2)
(
(*stcl)[1]=p;(*stcl)[2]=q;(*stcl)[3]=p+1;(*stcl)[4]=q;(*stcl)[5]=p+2;(*stcl)[6]=q;
(*stcl)[7]=p-1;(*stcl)[8]=q; (*stcl)[9]=p-2;(*stcl)[10]=q;(*stcl)[11]=p;(*stcl)[12]=q-1;
(*stcl)[13]=p;(*stcl)[14]=q-2;(*stcl)[15]=p;(*stcl)[16]=q+1;(*stcl)[17]=p;(*stcl)[18]=q+2;
(*stcl)[19]=p+1;(*stcl)[20]=q+1;(*stcl)[21]=p+1;(*stcl)[22]=q-1;(*stcl)[23]=p-1;
(*stcl)[24]=q-1;(*stcl)[25]=p-1; (*stcl)[26]=q+1;
)
)
else if (p>(N-2) && q>2 && q<=M)
(
if(p==N-1 && q>2 && q<=(M-2))
(
(*stcl)[1]=p;(*stcl)[2]=q;(*stcl)[3]=p+1;(*stcl)[4]=q;(*stcl)[5]=1;(*stcl)[6]=q;
(*stcl)[7]=p-1;(*stcl)[8]=q; (*stcl)[9]=p-2;(*stcl)[10]=q;(*stcl)[11]=p;(*stcl)[12]=q-1;
(*stcl)[13]=p;(*stcl)[14]=q-2;(*stcl)[15]=p;(*stcl)[16]=q+1;(*stcl)[17]=p;(*stcl)[18]=q+2;
(*stcl)[19]=p+1;(*stcl)[20]=q+1;(*stcl)[21]=p+1;(*stcl)[22]=q-1;(*stcl)[23]=p-1;
(*stcl)[24]=q-1; (*stcl)[25]=p-1; (*stcl)[26]=q+1;
)
else if(p==N && q>2 && q<=(M-2))
(
(*stcl)[1]=p;(*stcl)[2]=q;(*stcl)[3]=1;(*stcl)[4]=q;(*stcl)[5]=2;(*stcl)[6]=q;
(*stcl)[7]=p-1;(*stcl)[8]=q; (*stcl)[9]=p-2;(*stcl)[10]=q;(*stcl)[11]=p;(*stcl)[12]=q-1;
(*stcl)[13]=p;(*stcl)[14]=q-2;(*stcl)[15]=p;(*stcl)[16]=q+1;(*stcl)[17]=p;
(*stcl)[18]=q+2;(*stcl)[19]=1;(*stcl)[20]=q+1;(*stcl)[21]=1;(*stcl)[22]=q-1;
(*stcl)[23]=p-1; (*stcl)[24]=q-1; (*stcl)[25]=p-1; (*stcl)[26]=q+1;
)
else if(p==(N-1) && q==(M-1))
(
(*stcl)[1]=p;(*stcl)[2]=q;(*stcl)[3]=p+1;(*stcl)[4]=q;(*stcl)[5]=1;(*stcl)[6]=q;
(*stcl)[7]=p-1;(*stcl)[8]=q; (*stcl)[9]=p-2;(*stcl)[10]=q;(*stcl)[11]=p;(*stcl)[12]=q-1;
(*stcl)[13]=p;(*stcl)[14]=q-2;(*stcl)[15]=p;(*stcl)[16]=q+1;(*stcl)[17]=-1;(*stcl)[18]=-1;
(*stcl)[19]=p+1;(*stcl)[20]=q+1;(*stcl)[21]=p+1;(*stcl)[22]=q-1;(*stcl)[23]=p-1;
(*stcl)[24]=q-1; (*stcl)[25]=p-1; (*stcl)[26]=q+1;
(
else if(p==N && q==(M-1))
(
(*stcl)[1]=p;(*stcl)[2]=q;(*stcl)[3]=1;(*stcl)[4]=q;(*stcl)[5]=2;(*stcl)[6]=q;
67
(*stcl)[7]=p-1;(*stcl)[8]=q; (*stcl)[9]=p-2;(*stcl)[10]=q;(*stcl)[11]=p;(*stcl)[12]=q-1;
(*stcl)[13]=p;(*stcl)[14]=q-2;(*stcl)[15]=p;(*stcl)[16]=q+1;(*stcl)[17]=-1;(*stcl)[18]=-1;
(*stcl)[19]=1;(*stcl)[20]=q+1;(*stcl)[21]=1;(*stcl)[22]=q-1;(*stcl)[23]=p-1;
(*stcl)[24]=q-1;(*stcl)[25]=p-1; (*stcl)[26]=q+1;
)
)
for(I=0;I<= M;I++)
(
if(I==(q-1))
(
printf(?r1=lf/t dta=lf/n?,r1,dta);
(*val)[1]=(6)/(pow(x,4))-(2)/((r1+I*x)*(r1+I*x)*pow(x,2))
+(6)/((r1+I*x)*(r1+I*x)*(r1+I*x)*(r1+I*x)*pow(dta,4))
+(8)/((r1+I*x)*(r1+I*x)*pow(x ,2)*pow(dta,2))
-(8)/((r1+I*x)*(r1+I*x)*(r1+I*x)*(r1+I*x)*pow(dta,2));
printf(?val[1]=lf/n?,(*val)[1]);
(*val)[2]=-(4)/(pow(x,4))-(4)/((r1+I*x)*pow(x,3))
+(1)/(x*(r1+I*x)*(r1+I*x)*(r1+I*x))-(4)/(x*(r1+I*x)
(r1+I*x)*(r1+I*x)*pow(dta,2))-(4)/((r1+I*x)*(r1+I*x)*pow(x,2)*pow(dta,2))
+(1)/((r1+I*x)*(r1+I*x)*pow(x,2));
printf(?val[2]=lf/n?,(*val)[2]);
(*val)[3]=(1)/(pow(x,4))+(2)/((r1+I*x)*pow(x,3));
printf(?val[3]=lf/n?,(*val)[3]);
(*val)[4]=-(4)/(pow(x,4))+(4)/((r1+I*x)*pow(x,3))
-1/(x*(r1+I*x)*(r1+I*x)*(r1+I*x))+(4)/(x*(r1+I*x)*
(r1+I*x)*(r1+I*x)*pow(dta,2))-(4)/((r1+I*x)*(r1+I*x)*pow(x,2)*pow(dta,2))
+1/((r1+I*x)*(r1+I*x)*pow(x,2));
printf(?val[4]=lf/n?,(*val)[4]);
(*val)[5]=1/(pow(x,4))+(2)/((r1+I*x)*pow(x,3));
printf(?val[5]=lf/n?,(*val)[5]);
(*val)[6]=-(4)/((r1+I*x)*(r1+I*x)*(r1+I*x)*(r1+I*x)*pow(dta,4))
-(4)/((r1+I*x)*(r1+I*x)*pow(x,2)*pow(dta,2))
+(4)/((r1+I*x)*(r1+I*x)*(r1+I*x)*(r1+I*x)*pow(dta,2))-(c*omega)/(2*dta);
printf(?val[6]=lf/n?,(*val)[6]);
(*val)[7]= 1/((r1+I*x)*(r1+I*x)*(r1+I*x)*(r1+I*x)*pow(dta,4));
printf(?val[7]=lf/n?,(*val)[7]);
(*val)[8]=-(4)/((r1+I*x)*(r1+I*x)*(r1+I*x)*(r1+I*x)*pow(dta,4))
-(4)/((r1+I*x)*(r1+I*x/2)*pow(x,2)*pow(dta,2))+(c*omega)/(2*dta);
printf(?val[8]= lf/n?,(*val)[8]);
(*val)[9]=(1)/((r1+I*x)*(r1+I*x)*(r1+I*x)*(r1+I*x)*pow(dta,4));
printf(?val[9]=lf/n?,(*val)[9]); (*val)[10]=(2)/(x*(r1+I*x)*(r1+I*x)*(r1+I*x)*pow(dta,2))
+(2)/((r1+I*x)*(r1+I*x)*pow(x,2)*pow(dta,2));
68
printf(?val[10]= lf/n?,(*val)[10]);
(*val)[11]=(2)/(x*(r1+I*x)*(r1+I*x)*(r1+I*x)*pow(dta,2))
+(2)/((r1+I*x)*(r1+I*x)*pow(x,2)*pow(dta,2));
printf(?val[11]=lf/n?,(*val)[11]); (*val)[12]=-(2)/(x*(r1+I*x)*(r1+I*x)*(r1+I*x)*pow(dta,2.0))
+(2)/((r1+I*x)*(r1+I*x)*pow(x,2.0)*pow(dta,2.0));
printf(?val[12]= lf/n?,(*val)[12]);
(*val)[13]=-(2)/(x*(r1+I*x)*(r1+I*x)*(r1+I*x)*pow(dta,2))
+(2)/((r1+I*x)*(r1+I*x)*pow(x,2)*pow(dta,2));
printf(?val[13]= lf/n?,(*val)[13]);
(*val)[14]=(rho*x)/(pow(deltaT,2));
printf(?val[14]= lf/n?,(*val)[14]);
)
)
(*val)[15]=(2*pow(dta,2)*pow(r2,2)+2*mu*pow(x,2))/(pow(x,2));
printf(?val[15]=lf/n?,(*val)[15]);
(*val)[16]=-(r2*pow(dta,2)*(r2+mu*x))/(pow(x,2));
printf(?val[16]=lf/n?,(*val)[16]);
(*val)[17]=-(r2*pow(dta,2)*(r2-mu*x))/(pow(x,2));
printf(?val[17]=lf/n?,(*val)[17]);
(*val)[18]=-1.00;
printf(?val[18]=lf/n?,(*val)[18]);
(*val)[19]=(-1/pow(x,2)-(mu/pow(r2,2.0)*x)+(2/pow(x,3))+(1/r2*pow(x,2))
-(1/pow(r2,2.0))-(2*(1-mu)/r2*x*pow(dta,2.0)))/(-2/pow(x,2.0)
-2*mu/pow(r2,2)*pow(dta,2.0)+2/r2*x-2*(1-mu)/pow(r2,3.0)*pow(dta,2.0));
printf(?val[19]=lf/n?,(*val)[19]);
(*val)[20]=(-1/pow(x,2.0)+(mu/pow(r2,2.0)*x)-(2/pow(x,3.0))
+(1/r2*pow(x,2.0))+(1/pow(r2,2.0))+(2*(1-mu)/r2*x*pow(dta,2.0)));
printf(?val[20]=lf/n?,(*val)[20]);
(*val)[21]=(-mu/pow(r2,2.0)*pow(dta,2.0)+(1/pow(r2,2)*pow(dta,2.0))
-((1-mu)/pow(r2,3.0)));
printf(?val[21]=lf/n?,(*val)[21]);
(*val)[22]=(-mu/pow(r2,2.0)*pow(dta,2.0)+(1/pow(r2,2)*pow(dta,2.0))
-((1-mu)/pow(r2,3.0)));
printf(?val[22]=lf/n?,(*val)[22]);
(*val)[23]=(1/pow(x,3.0));
printf(?val[23]=lf/n?,(*val)[23]);
(*val)[24]=(-1/pow(x,3.0));
printf(?val[24]=lf/n?,(*val)[24]);
(*val)[25]=((1-mu)/r2*x*pow(dta,2.0));
printf(?val[25]=lf/n?,(*val)[25]);
(*val)[26]=((1-mu)/r2*x*pow(dta,2.0));
printf(?val[26]=lf/n?,(*val)[26]);
69
(*val)[27]=-((1-mu)/r2*x*pow(dta,2.0));
printf(?val[27]=lf/n?,(*val)[27]);
(*val)[28]=-((1-mu)/r2*x*pow(dta,2.0));
printf(?val[28]=lf/n?,(*val)[28]);
)
70
Appendix B
Matlab program to predict the response of the disc
Matlab program to get the displacement in transverse direction clear all;
clc;
int i;
int a;
int b;
i=1;
c=0.0002*eye(25,25)
R=0.045:0.0036666666666666666667:.133;
THETA=0:2*pi/24:2*pi;
[r,theta]=meshgrid(R,THETA);
x=r.*cos(theta);
y=r.*sin(theta);
A=[0 0.0005 0 0 0 0 0.0005 0 0 0 0 0.0005 0 0 0 0 0 0 0 0 0 0.0005 0 0 0];
F=diag(A);
load new55.txt;
f=new55;
fprintf(?f?,f);
s=f-c;
z=inv(s)*F;
fprintf(?the value of x=f/ty=/f/tz=/f/n?,x,y,z);
surf(x,y,z);
grid on
figure
plot(z)
71
Appendix C
C-program to solve Coupled Differential Equation
include < stdio.h >
include < math.h >
double mass(double u,double v,double y1, double y2);
double full(double u,double v,double y1, double y2);
int main(void)
( FILE *sucfile11,*sucfile12,*sucfile13,*sucfile14;
doubley1,y2,u,v,y1new,y2new,h,m1,m2,m3,m4,p1,p2,p3,p4,k1,k2,k3,k4,l1,l2,l3,l4,unew,vnew,f,g;
int t=0;
sucfile11=fopen(?c:/success11/success1.txt?,?w?);
sucfile12=fopen(?c:/success11/success2.txt?,?w?);
sucfile13=fopen(?c:/success11/success3.txt?,?w?);
sucfile14=fopen(?c:/success11/success4.txt?,?w?);
y1=0;
y2=0;
u=0;
v=0;
for(t=0;t< 800;t++)
(
h=0.1;
m1=h*mass(u,v,y1,y2);
printf(?m1=lf/n?,m1);
p1=h*full(u,v,y1,y2);
printf (?check in main : d/t mass1=lf/t ?,t,mass(u,v,y1,y2));
k1=h*u;
l1=h*v;
m2=h*mass(u+m1/2.0,v+p1/2.0,y1+k1/2.0,y2+l1/2.0);
p2=h*full(u+m1/2.0,v+p1/2.0,y1+k1/2.0,y2+l1/2.0);
printf(?checkinmain: d/tmass2=lf/t?,t,mass(u+m1/2.0,v+p1/2.0,y1+k1/2.0,y2+l1/2.0));
k2=h*(u+m1/2);
l2=h*(v+p1/2);
m3=h*mass(u+m2/2.0,v+p2/2.0,y1+k2/2.0,y2+l2/2.0);
p3=h*full(u+m2/2.0,v+p2/2.0,y1+k2/2.0,y2+l2/2.0); k3=h*(u+m2/2.0);
l3=h*(v+p2/2.0);
m4=h*mass(u+m3,v+p3,y1+k3,y2+l3);
p4=h*full(u+m3,v+p3,y1+k3,y2+l3);
k4=h*(u+m3);
72
l4=h*(v+p3);
y1new=y1+k1*(1/6.0)+k2*(1/3.0)+k3*(1/3.0)+k4*(1/6.0);
y2new=y2+l1*(1/6.0)+l2*(1/3.0)+l3*(1/3.0)+l4*(1/6.0);
printf(?y1=lf /n?,y1new);
fprintf(sucfile11,?lf /n?,y1new);
fprintf(sucfile12,?lf /n?,y2new);
unew=u+m1*(1/6.0)+m2*(1/3.0)+m3*(1/3.0)+m4*(1/6.0);
vnew=v+p1*(1/6.0)+p2*(1/3.0)+p3*(1/3.0)+p4*(1/6.0);
fprintf(sucfile13,?lf/n?,unew);
fprintf(sucfile14,?lf/n?,vnew);
printf (?new values : d/tlf/tlf/tlf/tlf/n?,t,y1new, y2new, unew, vnew);
y1=y1new;
y2=y2new;
u=unew;
v=vnew;
)
fclose(sucfile11);
fclose(sucfile12);
fclose(sucfile13);
fclose(sucfile14);
)
double mass(double u,double v,double y1, double y2)
(
double Nstatic=10;
double M1=2.5;
double c=0.1;
double k=20;
double f;
f=((Nstatic-c*(u-v)-k*(y1-y2)-k*pow((y1-y2),3.0))/M1);
printf (?check in function mass= : lf/n?,f );
return f;
)
double full(double u,double v,double y1, double y2)
(
double M2=0.25;
double c=0.1;
double k=20;
double g;
g=((c*(u-v)+k*(y1-y2)+k*pow((y1-y2),3.0))/M2);
return g;
)
73
C-program to solve Coupled Differential Equation Continued
include < stdio.h >
include < math.h >
define N 10
double dev(double x1,double x2,double r,double s,double y1,double y2,double u,
double v);
double devi(double x1,double x2,double r,double s,double y1,double y2,double u,
double v);
int main()
(
FILE *successfile11,*successfile12,*successfile13,*successfile14,
gsuccessfile11,*gsuccessfile12,*gsuccessfile13,*gsuccessfile14;
int t=0;
double y1[800],y2[800],u[800],v[800];
gsuccessfile11=fopen(?c://success11//success01.txt?,?w?);
gsuccessfile12=fopen(?c://success11//success02.txt?,?w?);
gsuccessfile13=fopen(?c://success11//success03.txt?,?w?);
gsuccessfile14=fopen(?c://success11//success04.txt?,?w?);
successfile11=fopen(?c://success11//success1.txt?,?r?);
successfile12=fopen(?c://success11//success2.txt?,?r?);
successfile13=fopen(?c://success11//success3.txt?,?r?);
successfile14=fopen(?c://success11//success4.txt?,?r?);
for (t=0;t< 800;t++)
(
fscanf(successfile11,?lf?,y1[t]);
fscanf(successfile12,?lf?,y2[t]);
fscanf(successfile13,?lf?,u[t]);
fscanf(successfile14,?lf?,v[t]);
)
fclose(successfile11);
fclose(successfile12);
fclose(successfile13);
fclose(successfile14);
)
double h,k1[800],k2[800],k3[800],k4[800],l1[800],l2[800],l3[800],
l4[800],x1[800],x2[800],r[800],s[800];
double m1[800],m2[800],m3[800],m4[800],p1[800],p2[800],p3[800],p4[800];
r[0]=0;
s[0]=0;
x1[0]=0;
x2[0]=0;
for(t=0;t< 800;t++)
74
(
h=0.1;
printf(?recheck the values:lf/nlf/n?,y1[t],y2[t]);
k1[t]=h*r[t];
printf(?k1=lf/n?,k1[t]);
l1[t]=h*s[t];
printf(?l1=lf/n?,l1[t]);
m1[t]=h*dev(x1[t],x2[t],r[t],s[t],y1[t],y2[t],u[t],v[t]);
printf(?m1=lf/n?,m1[t]);
p1[t]=h*devi(x1[t],x2[t],r[t],s[t],y1[t],y2[t],u[t],v[t]);
printf (?check in main d/tdev1=lf/tdevi1=lf/n?,t,dev(x1[t],x2[t],r[t],s[t],
y1[t],y2[t],u[t],v[t]),devi(x1[t],x2[t],r[t],s[t],y1[t],y2[t],u[t],v[t]));
k2[t]=h*(r[t]+0.5*m1[t]);
l2[t]=h*(s[t]+0.5*p1[t]);
m2[t]=h*dev(x1[t]+0.5*k1[t],x2[t],r[t]+0.5*m1[t],s[t],y1[t],y2[t],u[t],v[t]);
p2[t]=h*devi(x1[t],x2[t]+0.5*l1[t],r[t],s[t]+0.5*p1[t],y1[t],y2[t],u[t],v[t]);
k3[t]=h*(r[t]+0.5*m2[t]);
l3[t]=h*(s[t]+0.5*p2[t]);
m3[t]=h*dev(x1[t]+0.5*k2[t],x2[t],r[t]+0.5*m2[t],s[t],y1[t],y2[t],u[t],v[t]);
p3[t]=h*devi(x1[t],x2[t]+0.5*k2[t],r[t],s[t]+0.5*p2[t],y1[t],y2[t],u[t],v[t]);
k4[t]=h*(r[t]+m3[t]);
l4[t]=h*(s[t]+p3[t]);
m4[t]=h*dev(x1[t]+ k3[t],x2[t],r[t]+ m3[t],s[t],y1[t],y2[t],u[t],v[t]);
p4[t]=h*devi(x1[t],x2[t]+ k3[t],r[t],s[t]+ p3[t],y1[t],y2[t],u[t],v[t]);
printf(?check: lf/nlf/n?,k1[t],k2[t]);
x1[t+1]=(x1[t]+0.167*k1[t]+0.334*k2[t]+0.334*k3[t]+0.167*k4[t]);
x2[t+1]=(x2[t]+0.167*l1[t]+0.334*l2[t]+0.334*l3[t]+0.167*l4[t]);
s[t+1]=s[t]+0.167*p1[t]+0.334*p2[t]+0.334*p3[t]+0.167*p4[t];
printf(?check the values lf/nlf/nlf/nlf/n?,x1[t+1],x2[t+1],r[t+1],s[t+1]);
fprintf(gsuccessfile11,?lf/n?,x1[t+1]);
fprintf(gsuccessfile12,?lf/n?,x2[t+1]);
fprintf(gsuccessfile13,?lf/n?,r[t+1]);
fprintf(gsuccessfile14,?lf/n?,s[t+1]);
x1[t]=x1[t+1];
x2[t]=x2[t+1];
r[t]=r[t+1];
s[t]=s[t+1];
)
fclose(gsuccessfile11);
fclose(gsuccessfile12);
fclose(gsuccessfile13);
fclose(gsuccessfile14);
75
)
double dev(double x1,double x2 ,double r,double s,double y1 ,double y2,double u,
double v)
(
double Nstatic=;
double M1=2.5;
int c1=1,c=1,c2=0.01;
double k1=0.01,k=20,k2=1;
double w,vr,Fric1,Fric2,Ff;
double vo=1;
double ud,us=0.3,b=0.12;
vr=vo+s-r;
ud=us-b*vr;
scanf(?vo=lf/n?,vo);
vr=vo+s-r;
printf(?vr=lf/n?,vr);
Fric1=-k1*x1-c1*r+k2*x2+c2*x2;
if (Fric1 < 0)
(
Fric1= -Fric1;
)
else
(
Fric1= Fric1;
)
printf(?Fric1.1=lf/n?,Fric1);
Fric2=us*Nstatic;
if (Fric2 < 0)
(
Fric2= -Fric2;
)
else
(
Fric2= Fric2;
)
if (vr==0)
(
if (Fric1 ? Fric2 && Fric1 ? 0)
(
Ff=-Fric1;
)
else if (Fric1 < Fric2 && Fric1 ? 0)
76
(
Ff=Fric1;
)
else if (Fric1 > Fric2 && Fric2 ? 0)
(
Ff=Fric2;
)
else if (Fric1 > Fric2 && Fric2 ? 0)
(
Ff=-Fric2;
)
w=(Ff-c1*r-k1*x1)/M1;
)
else if(vr > 0)
(
w=(-ud*(Nstatic+c*(u-v)+k*(y1-y2)+k*pow((y1-y2),3.0))-c1*r-k1*x1)/M1;
printf(?slip001=lf/n?,w);
)
else
(
w=(ud*(Nstatic+c*(u-v)+k*(y1-y2)+k*pow((y1-y2),3.0))-c1*r-k1*x1)/M1;
printf(?slip002=lf/n?,w);
)
return w;
)
double devi(double x1,double x2,double r, double s, double y1, double y2,double u,
double v)
(
double Nstatic=10;
double M2=1;
double c1=.01,c=0.1,c2=0.01;
double k1=1,k=20,k2=1;
double z,vr,Fric1,Fric2;
double vo=1;
double us=0.3;
double ud,Ff,b=0.12;
vr=vo+s-r;
ud=us-b*vr;
Fric1=-k1*x1-c1*r+k2*x2+c2*x2;
printf(?Fric2=lf/n?,Fric1);
if (Fric1 < 0)
77
(
Fric1= -Fric1;
)
else
(
Fric1= Fric1;
)
printf(?Fric1=lf/n?,Fric1);
Fric2=us*Nstatic;
if (Fric2 < 0)
(
Fric2= -Fric2;
)
else
(
Fric2= Fric2;
)
printf(?Fric2.2=lf/n?,Fric2);
if (vr==0)
(
if (Fric1 ? Fric2 && Fric1 ? 0)
(
Ff=-Fric1;
)
else if (Fric1 < Fric2 && Fric1 ? 0)
(
Ff=Fric1;
)
else if (Fric1 > Fric2 && Fric2?0)
(
Ff=Fric2;
)
else if (Fric1 > Fric2 && Fric2 ? 0)
(
Ff=-Fric2;
)
z=(-Ff-c2*s-k2*x2)/M2;
printf(?stick01=lf/n?,z);
)
else if (vr > 0)
(
z=(-ud*(Nstatic+c*(u-v)+k*(y1-y2)+k*pow((y1-y2),3.0)-c2*s-k2*x2))/M2;
78
printf(?slip01=lf/n?,z);
)
else
(
z=ud*(Nstatic+ c*(u-v)+ k*(y1-y2)+ k*(y1-y2)*(y1-y2)*(y1-y2))/M2;
printf(?slip02=lf/n?,z);
)
return z;
)
79