Modelling, Simulation and Validation of Dynamic

Systems with Lab-Experiments

Helmut Scherf

Exam: test 90 min without expedients

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-1-
 Modelling and Model: Some Statements

Compact summary of knowledge

Easier or safer or less expensive to work with than the real world
Sometimes a necessity because the system does not exist
Modelling for control
Modelling for design
Modelling for education
A model only reflects a (relevant) selection of the original’s properties
Try to assess the validity and the accuracy of the model
– Assumptions made, effects that are neglected
– Accuracy and uncertainty
Einstein: A model should be as simple as possible but not simpler!
Ralph P. Schlenker, Exxon, 1998:
Modelling and simulation technologies are key to achieving
manufacturing excellence and for assessing risk in unit operations.
As we make our plant more flexible to respond to business
opportunities efficient modelling and simulation techniques will
become commonly used tools.

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-2-
 Modelling and Model

S. W. Golomb, Simulation 14 (1970), 197-198 : Some Don’ts!

DON‘T believe that the model is the reality

DON‘T extrapolate beyond the region of fit

DON‘T distort reality to fit the model

DON‘T retain a discredited model

DON‘T fall in love with your model

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-3-
 Modelling and Simulation of Dynamic Systems

System dynamics is dealing with the value of states in a system over time.
Set up of differential equation(s) of a dynamic system
Static systems are described by algebraic equations

•Mass balance
•Flow and heat balance
•laws of thermodynamics,
•Principle of energy conservation
•Newton’ s law Model (white box)
•Principle of linear momentum Opposite: black box
•Lagrange equations
•Kirchhoff’s voltage law
•Kirchhoff’s current law
•etc.

Simulation with MATLAB &Simulink

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-4-
 Simulation

(Shannon, 1975)

Simulation is the process of designing a model of a real

system and conducting experiments with this model for
the purpose either of understanding the behavior of the
system and its underlying causes or of evaluating
various designs of an artificial system or strategies for
the operation of the system.

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-5-
 Simulation Process
Problem

Set up Differential
Equations

Determine System
Parameters

Simulation

Comparison
Simulation no Model
Measurement adaptation
ok?
yes
Stop
Fakultät für Maschinenbau und Mechatronik
Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-6-
 Oscillator: Real Model

Copper plate (mass)

Leaf spring
Ultra sonic
sensor

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-7-
 Oscillator: Equivalent Model

coil spring
mass
k

m
d

damper
frictionless

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-8-
 Oscillator: Free Body Diagram
Spring relaxed x
Static position of rest
k

m
d

mass deflected
k x spring extended

FS Fi
m
d FD

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-9-
 Oscillator: Set-up of the Differential Equation

spring force FS = k ⋅ x
damping force FD = d ⋅ x&
inertial force Fi = m ⋅ &x&

Equilibrium of forces: Newton’s law or d’Alembert’s principle

Differential equation: m ⋅ &x& + d ⋅ x& + k ⋅ x = 0

Initial conditions: x(0) = x0 , x& (0) = 0

•ordinary
•homogeneous (free oscillation) <-> forced oscillation, inhomogeneous diff. eq.
•second order (2 energy storages)

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-10-
 Differential Equation

m d
m ⋅ &x& + d ⋅ x& + k ⋅ x = 0 → ⋅ &x& + ⋅ x& + x = 0
k k
Alternative form: T 2 ⋅ &x& + 2 ⋅ D ⋅ T ⋅ x& + x = 0
Equating the corresponding coefficients results:

m d
T = , 2 ⋅ D ⋅T =
2

k k
1
T= time constant [sec]
ω0
2π
ω0 = angular frequency of the undamped system [rad/sec]
TC TC period (cycle duration) [sec]
D dimensionless damping coefficient [-]

ωd = ω0 1 − D 2 angular frequency of the underdamped system [rad/sec]

Fakultät für Maschinenbau und Mechatronik
Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-11-
 Solution of the Differential Equation

Exponential approach: x(t ) = e λ ⋅t

Second derivative: &x& = λ2 ⋅ e λ ⋅t First derivative: x& = λ ⋅ e λ ⋅t

Plug in: T 2 ⋅ &x& + 2 ⋅ D ⋅ T ⋅ x& + x = 0 →

T 2 ⋅ λ2 ⋅ e λ ⋅t + 2 ⋅ D ⋅ T ⋅ λ ⋅ e λ ⋅t + e λ ⋅t = 0
2D 1
Characteristic equation: ⇒λ + ⋅λ + 2 = 0
2

T T
quadratic equation ⇒ two solutions

D D2 1 D D2 −1
λ1/ 2 =− ± 2
− 2 =− ±
T T T T T
Fakultät für Maschinenbau und Mechatronik
Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-12-
 Solution of the Differential Equation
Four cases: overdamped D >1
critically damped D =1
underdamped 0 < D <1
(undamped) D=0

If the roots are: General solution:

distinct x(t ) = C1 ⋅ e λ1 ⋅t + C2 ⋅ e λ2 ⋅t
Calculation of arbitrary constants
λ ⋅t λ ⋅t
repeated, single x(t ) = C1 ⋅ e + C2 ⋅ t ⋅ e C1 and C2 with initial conditions!
roots

Complex conjugate x(t ) = eα ⋅t ⋅ (C1 ⋅ cos β ⋅ t + C2 ⋅ sin β ⋅ t )

λ1/ 2 = α ± j ⋅ β

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-13-
 Solution of the Differential Equation

Here: underdamped system: 0 < D < 1⇒ λ complex

D D2 1 D 1− D2
λ1/ 2 =− ± − = − ± j ⋅ = − D ⋅ ω ± j ⋅ ω ⋅ 1 − D 2

T T2 T2 T T 1
424 30 10 4243
α β

General solution of the homogeneous differential equation:

x(t ) = eα ⋅t ⋅ (C1 ⋅ cos β ⋅ t + C2 ⋅ sin β ⋅ t )

(
= e − D⋅ω0 ⋅t ⋅ C1 ⋅ cos ω0 ⋅ 1 − D 2 ⋅ t + C2 ⋅ sin ω0 ⋅ 1 − D 2 ⋅ t )

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-14-
 Calculation of the coefficients

x(0) = x0 = C1

(
x& (t ) = − D ⋅ ω0 ⋅ e − D⋅ω0 ⋅t ⋅ C1 ⋅ cos ω0 ⋅ 1 − D 2 ⋅ t + C2 ⋅ sin ω0 ⋅ 1 − D 2 ⋅ t )
(
+ e − D⋅ω0 ⋅t ⋅ − C1 ⋅ ω0 ⋅ 1 − D 2 ⋅ sin ω0 ⋅ 1 − D 2 ⋅ t + C2 ⋅ ω0 ⋅ 1 − D 2 ⋅ cos ω0 ⋅ 1 − D 2 ⋅ t )
x& (0) = 0 = − D ⋅ ω0 ⋅ C1 + C2 ⋅ ω0 ⋅ 1 − D 2 = − D ⋅ x0 + C2 ⋅ 1 − D 2
D ⋅ x0
⇒ C2 =
1− D2
General solution of the homogeneous differential equation:

− D⋅ω0 ⋅t  D ⋅ x0 
x(t ) = e ⋅  x0 ⋅ cos ω0 ⋅ 1 − D ⋅ t +
2
⋅ sin ω0 ⋅ 1 − D ⋅ t 
2

 1− D 2


Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-15-
 Solving the Differential Equation with MATLAB
syms w0 d %don‘t use D, because D is the differential operator
S=dsolve('D2x/w0^2+2*d/w0*Dx+x=0')
d=2; %overdamped
S=dsolve('D2x/w0^2+2*2/w0*Dx+x=0')
d=1; %critically damped
S=dsolve('D2x/w0^2+2*1/w0*Dx+x=0')
d=0.6; %underdamped
S=dsolve('D2x/w0^2+2*0.6/w0*Dx+x=0')

%Now with initial conditions

S=dsolve('D2x/w0^2+2*0.6/w0*Dx+x=0','x(0)=x0','Dx(0)=0')
w0=0.2; x0=0.3;
S=dsolve('D2x/0.2^2+2*0.6/0.2*Dx+x=0','x(0)=0.3','Dx(0)=0')
subs(S) %replaces the symbolic variables by numeric values
%compare this solution with solution calculated by hand
9/40
0.6*x0/sqrt(1-0.6^2)
4/25
w0*sqrt(1-0.6^2)

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-16-
 Solving the Differential Equation with MATLAB
syms x(t) w0 d x0 %don‘t use D, because D is the differential operator
S=dsolve(diff(x,2)/w0^2+2*d/w0*diff(x)+x==0)
d=2; %overdamped
S=dsolve(diff(x,2)/w0^2+2*2/w0*diff(x)+x==0)
d=1; %critically damped
S=dsolve(diff(x,2)/w0^2+2*1/w0*diff(x)+x==0)
d=0.6; %underdamped
S=dsolve(diff(x,2)/w0^2+2*0.6/w0*diff(x)+x==0)
%Now with initial conditions
Dx=diff(x); D2x=diff(x,2)
S=dsolve(D2x/w0^2+2*0.6/w0*Dx+x==0, x(0)==x0, Dx(0)==0)
w0=0.2; x0=0.3;
S=dsolve(D2x/w0^2+2*0.6/w0*Dx+x==0, x(0)==x0, Dx(0)==0)
subs(S) %replaces the symbolic variables by numeric values
%compare this solution with solution calculated by hand
9/40
0.6*x0/sqrt(1-0.6^2)
4/25
w0*sqrt(1-0.6^2)

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-17-
 Types of Damping

Different types of damping:

•Viscous damping: damping force is proportional to the velocity

m ⋅ &x& + d ⋅ x& + k ⋅ x = 0, x(0) = x0 , x& (0) = 0

•Sliding friction damping: absolute value of damping force is constant, sign!

m ⋅ &x& + d f ⋅ sign( x& ) + k ⋅ x = 0, x(0) = x0 , x& (0) = 0

•quadratic damping: damping force is proportional to the velocity squared,

sign!

m ⋅ &x& + d q ⋅ sign( x& ) ⋅ x& 2 + k ⋅ x = 0, x(0) = x0 , x& (0) = 0

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-18-
 Types of Damping
viscous damping
1 −d
⋅t
2⋅m Exponential decrease
e
0

-1
0 5 10 15 20 25 30 35 40 45 50
time in sec
friction damping
1

Linear decrease
0

-1
0 5 10 15 20 25 30 35 40 45 50
time in sec
air damping
1
Strong decrease weak decrease
0

-1
0 5 10 15 20 25 30 35 40 45 50
time in sec

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-19-
 Block Diagram

viscous (magnetic) damping (linear dependence)

x'' x' x
1 1
1
s s

1/mass Integrator Integrator 1 Scope

d damping coefficient

0.1

k spring constant

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-20-
 Block Diagram

friction damping
x '' x' x
1 1
1
s s

Integrator 2 Integrator 3 Scope

1 / mass 1
d damping coefficient 1 Sign

0 . 05

k spring constant 1

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-21-
 Block Diagram

quadratic (air ) damping (quadratic dependence)

x'' x' x
1 1
1
s s

Integrator 4 Integrator 5 Scope

1 / mass 2
Fcn

sgn (u )* 0 . 4 * u ^2

k spring constant 2

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-22-
 How to determine the System Parameters

Measuring the mass m with a balance m = 0.932 kg

Determination of the spring constant k : 2 methods

First method: rotate through 90° and measure the displacement of the mass

x = 11 mm

F
k = = 831 N/m
x

F = m⋅ g
Fakultät für Maschinenbau und Mechatronik
Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-23-
 How to determine the System Parameters

Second method: beam theory with concentrated force on free end

(strength theory)

cantilever Bending beam

h F deflection
beam with F ⋅l3
fixed end x = 4⋅
E ⋅ b ⋅ h3
l
E-modulus
b
width b = 7 mm E ⋅ b ⋅ h3
⇒ k = 2⋅
length l = 88 mm 4⋅l3
thickness h = 0.8 mm

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-24-
 How to determine the System Parameters
 xˆ1  1  xˆ1 
ϑ = ln  = ln 
Determination of the damping constant d :  2  n  xˆn 
ˆ
x
logarithmic damping decrement ϑ ϑ
D=
viscous damping 4 ⋅π 2 + ϑ 2
1

x̂1 x̂n
x̂2
0

-1
0 5 10 15 20 25 30 35 40 45 50

time in sec

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-25-
 Block Diagram

magnetic damping

x'' x' x
1 1
1/m s s

Gain 2 Integrator Integrator 1 Scope

Gain

Gain 1

m ⋅ &x& + d ⋅ x& + k ⋅ x = 0, x(0) = x0 , x& (0) = 0

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-26-
 MATLAB Code
%Oscillator with magnetic damping
%--------------------------------
load [Link]
t=MagneticDamping(:,1); %time vector [s]
xm=MagneticDamping(:,2); %measured voltage ultrasonic sensor [V]
plot(t-0.9,xm-2.75)
grid
axis([0 6 -1 1])
theta=log(0.776/0.67); %logarithmic damping decrement [-]
D=theta/sqrt(4*pi^2+theta^2) %dimensionless damping coefficient [-]

m=0.932; %mass [kg]

g=9.81;
k=m*g/11e-3; %spring constant [N/m] See remark!
b_spring=7e-3; %width [m]
l_spring=88e-3; %length [m]
h_spring=0.8e-3; %thickness [m]
E=2.7e+011; %E-modulus spring steel [N/m^2]
k=E*b_spring*h_spring^3/2/l_spring^3 % See remark!
See remark!

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-27-
 MATLAB Code
x0=-0.97; %initial value displacement integrator [V]
w0=sqrt(k/m); %angular frequency of the undamped system [rad/sec]
d=2*D*k/w0; %viscous damping coefficient [N/(m/s)]
wd=w0*sqrt(1-D^2); %angular frequency of the underdamped system [rad/sec]

TC=2*pi/wd; %cycle duration [sec]

sim('oscillatormagneticdamping')
plot(t-0.9,xm-2.75,x(:,1),x(:,2),'r')
grid
axis([0 6 -1 1])

Remark:
different methods calculating the spring constant lead to different results

S=dsolve('D2x/27.6006^2+2*0.0234/27.6006*Dx+x=0','x(0)=-0.97','Dx(0)=0')
t=0:0.01:6;
y=eval(S);
plot(t,y)
Fakultät für Maschinenbau und Mechatronik
Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-28-
 Comparison of Simulation and Measurement

1
measurement
simulation
0.8

0.6

0.4

0.2
displacement

-0.2

-0.4

-0.6

-0.8

-1
0 1 2 3 4 5 6
time in sec

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-29-
 Block Diagram

friction damping
x'' 1 x' 1 x
1/m
s s

1/mass Integrator 2 Integrator 3 Scope

d damping coefficient Sign

26

k spring constant

-K-

m ⋅ &x& + d f ⋅ sign( x& ) + k ⋅ x = 0, x(0) = x0 , x& (0) = 0

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-30-
 MATLAB Code
%Oscillator with friction damping
%--------------------------------
load [Link]
t=FrictionDamping(:,1); %time vector [s]
xm=FrictionDamping(:,2); %measured voltage ultrasonic sensor [V]
plot(t-0.62,xm-2.82)
grid
axis([0 3.5 -2 2])
m=0.932; %mass [kg]
%take over the values
b_spring=7e-3; %width [m]
l_spring=88e-3; %length [m]
h_spring=0.8e-3; %thickness [m]
E=2.7e+011; %E-modulus spring steel [N/m^2]
k=E*b_spring*h_spring^3/2/l_spring^3 %spring constant [N/m]
x0=-1.68; %initial value displacement integrator [V]
sim('oscillatorfrictiondamping')
plot(t-0.62,xm-2.82,x(:,1),x(:,2),'r')
grid
axis([0 3.5 -2 2])
xlabel('time in sec')
ylabel('displacement')
legend('measurement','simulation')

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-31-
 Comparison of Simulation and Measurement

2
measurement
simulation
1.5

0.5
displacement

-0.5

-1

-1.5

-2
0 0.5 1 1.5 2 2.5 3 3.5
time in sec

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-32-
 Block Diagram

quadratic (air) damping

x'' 1 x' 1 x
1/m
s s

1/mass2 Integrator 4 Integrator 5 Scope

Fcn

sgn(u)*0.015 *u^2

k spring constant 2

-K-

m ⋅ &x& + d q ⋅ sign( x& ) ⋅ x& 2 + k ⋅ x = 0, x(0) = x0 , x& (0) = 0

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-33-
 MATLAB Code
%Oscillator with quadratic (air) damping
%--------------------------------
load [Link]
t=AirDamping(:,1); %time vector [s]
xm=AirDamping(:,2); %measured voltage ultrasonic sensor [V]
plot(t-1.26,xm-2.8)
grid
axis([0 40 -1.5 1.5])
m=0.932; %mass [kg]
g=9.81; %gravitational constant [m/s^2]
%take over the values
b_spring=7e-3; %width [m]
l_spring=88e-3; %length [m]
h_spring=0.8e-3; %thickness [m]
E=2.7e+011; %E-modulus spring steel [N/m^2]
k=E*b_spring*h_spring^3/2/l_spring^3 %spring constant [N/m]
x0=-1.2; %initial value displacement integrator [V]
sim('oscillatorairdamping')
plot(t-1.26,xm-2.8,x(:,1),x(:,2),'r')
grid
axis([0 40 -1.5 1.5])
xlabel('time in sec')
ylabel('displacement')
legend('measurement','simulation')
Fakultät für Maschinenbau und Mechatronik
Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-34-
 Comparison of Simulation and Measurement

1.5
measurement
simulation

0.5
displacement

-0.5

-1

-1.5
0 5 10 15 20 25 30 35 40
time in sec

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-35-
 DC-Motor

Tacho
DC-Motor
generator

Banana jack

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-36-
The “Anatomy” of a Permanent Excitation DC-Motor

carbon armature
commutator armature winding permanent
brushes
magnet

N S

uM

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-37-
 Torque Generation of a DC-Motor
Force F on a current-carrying conductor

F = B ⋅l ⋅i
current
N length
Magnetic flux density

i •Torque M = F ⋅ r = kG ⋅ i
•uM motor voltage
•i armature current

brush
S

i commutator

uM
+ -
Fakultät für Maschinenbau und Mechatronik
Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-38-
 Equivalent Circuit Diagram of a DC-Motor

i
R R ⋅i voltage drop across
the resistor

motor voltage uM
di
L L⋅ voltage drop across
the inductor
dt

= eM = kG ⋅ ω
Induced voltage is
proportional to angular
speed

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-39-
 Free Body Diagram of the Rotor

ϕ J ⋅ ϕ&&
MR
ML
MM

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-40-
 Differential Equations of a DC-Motor

di
Kirchhoff‘s voltage law: u M = R ⋅ i + L ⋅ + kG ⋅ ω
dt

Newton‘s law:
equilibrium of moments : M M − M R − M L = J ⋅ ϕ&&

Motor torque: M M = kG ⋅ i
MM : motor torque [Nm] kG : electromotive force constant [Vs],
torque constant [Nm/A]
ML : load torque [Nm]

MR : Friction torque [Nm]

J: mass moment of inertia [kg m2]

ϕ&& : angular acceleration [rad/sec]

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-41-
 Block Diagram

R=4.64; %Armature resistance [Ohm]

Step MLoad KG=0.0159; %Electromotive force constant [Vs
MLoad in Nm J=8.5e-7; % moment of inertia [kgm^2]
uTacho in V MR=2.4e-3; % friction torque [Nm]
uM in V
Scope
KGT=1419; %tachogenerator constant [min^-1/V
Step uM DC-Motor
0->5 V

1
MLoad in Nm
angular speed
i in A M in Nm rad/sec
1
2 1/R KG 1/J 60/2/pi 1/KGT 1
s
uM in V uTacho in V
Integrator
Gain 4 Gain 1 Gain Gain 5 Gain 3

MR Constant

Gain 2

KG

File name: [Link]

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-42-
 Determination of the System Parameters
1. Electromotive force constant kG

dϕ
eM = kG ⋅
dt
t

t ϕ ∫e M dt
∫ eM dt = kG ⋅ ∫ dϕ ⇒ kG =
0 0
0
ϕ

Connect the banana jack of the DC-motor with the analog input
Turn the motor shaft e. g. 10 times 2 π t
Calculate the constant with above equation eM dt ∫
kG = 0

Analog 1 ϕ
Input s
0.999Vsec
Analog Input Integrator Display
= = 0.0159 Vsec
0.015
Offset 10 ⋅ 2 ⋅ π
correction

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-43-
 Determination of the System Parameters
…..or with an analog integrator:
Problem! Drift of the integrator! ==>remedy: drift compensation

CI

RI -
t
uin = eM +
uout =−
1
⋅ ∫ em dt
RI ⋅ C I 0

RI = 100 kΩ
C I = 1µF
Fakultät für Maschinenbau und Mechatronik
Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-44-
 Determination of the System Parameters
…..or acquire the induced voltage with a storage scope and evaluation with MATLAB

0.7

0.6

0.5

0.4
induced voltage in V

0.3

0.2

0.1

-0.1
-20 -15 -10 -5 0 5 10 15 20
time in sec

load [Link]
m=10; %number of turns
sum(kgMeasuring(1:2047,2).*diff(kgMeasuring(:,1)))/(2*pi*m)

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-45-
 Determination of the System Parameters
2. Electric resistance R

Let the DC-motor run with e. g. u M = 4.9 V

di
Measure the current i in the steady state ⇒ L⋅ = 0
Measure the angular speed ω in the steady state dt
Tachogenerator constant k = 1419 min-1/V
T

di
u M = R ⋅ i + L ⋅ + kG ⋅ ω
dt
1
ω = 2 ⋅ π ⋅ n = 2 ⋅ π ⋅ kT ⋅ uT ⋅ = 2 ⋅ π ⋅1419/60 ⋅1.77 V = 263 rad/sec
60
rad
4.9V − 0.0159Vsec ⋅ 263
u M − kG ⋅ ω sec = 4.64 Ω
⇒R= =
i 0.155A

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-46-
 Determination of the System Parameters
3. Friction torque MR

Let the DC-motor run with u M = 1 V , 2 V, 3 V, 4 V, 5 V

Measure the current i in the steady state ⇒ J ⋅ ϕ&& = 0 ⇒ M M = M R
Measure the motor speed n in the steady state M R = kG ⋅ i
uM [ V] uT [V] i [mA]
1.6 0 190
1.6 0.28 146
2 0.482 143
3 0.92 146
4 1.37 152
5 1.79 157

Conclusion: the friction torque is approximately constant (dynamic friction)

Fakultät für Maschinenbau und Mechatronik
Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-47-
 Friction Torque
Breakaway x 10
-3

torque
3

2.5

Average torque M R = 2.4 ⋅10 −3 Nm

Friction Torque [Nm]

1.5

0.5

0
0 500 1000 1500 2000 2500 3000

Motor Speed [min-1]

Conclusion: the friction torque is approximately constant (dynamic friction)
Fakultät für Maschinenbau und Mechatronik
Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-48-
 Determination of the System Parameters
4. Moment of inertia J

Let the DC-motor run with uM = 5 V

At t =0 disconnect the DC-motor

Measure the angular speed ω with a scope

3.
M M − M R = J ⋅ ϕ&&
− MR − MR
t = 0 : − M R = J ⋅ ϕ&& ⇒ J = ≈
ϕ&& ∆ω
Scope
∆t (s. following slide)

− M R − 2.4 ⋅10 −3 Nm
J≈ = = 8.5 ⋅10 −7 kgm 2
∆ω rad
− 267
∆t sec
94 ⋅10 −3 sec
Fakultät für Maschinenbau und Mechatronik
Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-49-
 Run Down Experiment

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-50-
 Measurement of the Step Response

2
1.8 V

1.5

0.63*1.8 V=1.14 V
1
uT [V]

0.5

TM = 14 ⋅10−3 sec
-0.5
-0.02 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18
time [sec]

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-51-
Step Response: Comparison of Simulation and Measurement

1.8

1.6 measurement
simulation
1.4

1.2
uT [V]

0.8

0.6

0.4

0.2

0
-0.02 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18
time [sec]

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-52-
 Time Domain, Laplace Domain

•Looking at the step response: we can approximate the DC-motor by a first order system

•Differential equation: T ⋅ u&T + uT = K ⋅ u M

•Time constant T = 14 m sec

1,8V
•Gain: K= = 0,36
5V
•Using Laplace Transformation, the above modelling equations can be expressed in terms of s.
T ⋅ s ⋅U T ( s) + U T ( s) = K ⋅U M ( s)

U T ( s) K
•Transfer function: G(s) = =
U M (s) T ⋅ s + 1

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-53-
 Measuring the Frequency Response

)
output signal Input signal u M = u M ⋅ sin(ω ⋅ t )
(green) (blue)
)
uT = uT ⋅ sin(ω ⋅ t + ϕ ) frequency

Phase angle amplitude

amplitude
Same frequency
Fakultät für Maschinenbau und Mechatronik
Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-54-
 Measuring the Frequency Response
1.5

tϕ
1
u M (t )

0.5
u T (t ) û M
û T
0

-0.5

-1
settling process TP
-1.5
0 1 2 3 4 5 6 7 8 9 10
Time in sec

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-55-
 Measuring the Frequency Response

•Move the system in the operation point, if necessary nonlinear system

•Excite the system with a sinusoidal signal, a linear system responds with a sinusoidal signal
with the same frequency
•Wait until the settling process is finished, i. e. decay of the homogeneous solution
•Measure the output amplitude and the input amplitude with a dual channel scope
•Calculate the ratio of output amplitude to input amplitude (magnitude)
yˆ
G =
uˆ
 yˆ 
G dB = 20 ⋅ log10  
 uˆ  ϕ t
•Measure the phase shift between output and input ϕ = − ⋅ 360° = −tϕ ⋅ f ⋅ 360°
TP

output lags the input, the angle becomes

negative

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-56-
 Screen Shot at a low Frequency

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-57-
 Screen Shot at a high Frequency

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-58-
 Frequency Response: Measured Data

ω [rad/sec] uT peak-peak [V] ϕ [msec]

tϕ
1 3.55 0
3 3.58 8
7 3.54 20
20 3.43 14
40 3.1 14
60 2.71 11.2
100 2.18 9.7
300 0.83 4.96

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-59-
 Calculating the Frequency Response with MATLAB
%Measuring the frequency response
%---------------------------------
w=[1 3 7 20 40 60 100 300]; %angular speed [rad/sec]
f=w/2/pi; %frequency [Hz]
uM=[10 10 10 10 10 10 10 10]; %motor voltage peak-peak [V]
uT=[3.55 3.58 3.54 3.43 3.1 2.71 2.18 0.83]; %tacho voltage
peak-peak [V]
tphi=[0 8 20 14 14 11.2 9.7 4.96]/1000; %phase displacement
[sec]

phi=-tphi.*f*360; %phase angle (°]

FdB=20*log10(uT./uM); %magnitude [dB]

KM=1.8/5; %system gain from step response

TM=14e-3; %measured time constant [sec]
G=KM./(TM*j*w+1); %first order system
GdB=20*log10(abs(G));
psi=angle(G)*180/pi;

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-60-
Calculating the Frequency Response with MATLAB

%Plot commands
%--------------
subplot(2,1,1)
semilogx(w,FdB,w,GdB,'r')
grid
xlabel('frequency [rad/sec]')
ylabel('magnitude [dB]')
subplot(2,1,2)
semilogx(w,phi,w,psi,'r')
grid
xlabel('frequency [rad/sec]')
ylabel('phase [°]')

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-61-
 Bode Diagram
-5

[dB]
-10
magnitude

-15

measurement
-20 simulation: first order system
-25
0 1 2 3
10 10 10 10
frequency [rad/sec]

-20

-40
phase [°]

-60

-80

-100
0 1 2 3
10 10 10 10
frequency [rad/sec]

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-62-
 Tachogenerator with Filter Capacitor
Tacho
DC-Motor
generator

Tachogenerator signal noisy due to commutator

==>Filter capacitor for smoothing the signal
C= 205µF

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-63-
 Equivalent Circuit Diagram

L R

uC
uT = C

Tachogenerator

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-64-
 Differential Equation Transfer Function

di
Kirchhoff‘s voltage law: uT = R ⋅ i + L ⋅ + uC
dt

di
i = C ⋅ u&C ⇒ = C ⋅ u&&C
dt

Differential equation: uT = R ⋅ C ⋅ u&C + L ⋅ C ⋅ u&&C + uC

UC 1
Transfer function:
G ( s) = =
UT L ⋅ C ⋅ s 2 + R ⋅ C ⋅ s + 1
1
≈ Additional low pass filter
R ⋅C ⋅ s +1
Fakultät für Maschinenbau und Mechatronik
Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-65-
 Equivalent Circuit Diagram

Voltage divider: complex

impedance
UC ZC
=
UT R + Z L + ZC
Z L = j ⋅ω ⋅ L
1
j ⋅ω ⋅ C
=
1
R + j ⋅ω ⋅ L +
= uT
R
j ⋅ω ⋅ C
U C ( s = jω ) 1
=
U T ( s = jω ) L ⋅ C ⋅ s 2 + R ⋅ C ⋅ s + 1
1
ZC =
j ⋅ω ⋅ C uC
Fakultät für Maschinenbau und Mechatronik
Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-66-
 Data Sheet Tachogenerator

Tachogenerator constant

inductance

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-67-
Step Response: Comparison of Simulation and Measurement

1.8

1.6 measurement
simulation
1.4

1.2
uC [V]

0.8

0.6

0.4

0.2

0
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45
time [sec]

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-68-
 Water Tank

Plexiglas tube

pressure A
sensor AV
(4 to 20 mA
0 to 100 mbar p Ball valve
500 Ω Shunt) V&out

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-69-
 Differential Equation Water Tank

Volume balance:
In words:
dV
V&in − V&out = The volume flow rate in – volume flow rate out =
dt The time rate of change of fluid volume inside the tube

= AT ⋅ h&
AT : Cross sectional area of the tube
AV : Cross sectional area of the ball valve
g : Acceleration due to gravity
h : Water level

V&in = f (u P ) volume flow rate in is a function of the pump voltage u P

pump characteristic curve

V&out = AV ⋅ 2 ⋅ g ⋅ h
1424 3 volume flow rate out
Torricelli

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-70-
 Solving the Differential Equation
− AV ⋅ 2 ⋅ g ⋅ h = AT ⋅ h& Differential equation

− AV ⋅ 2 ⋅ g
k ⋅ h = h with k =
&
AT
1
−
k ⋅ dt = h dh 2 Separation of variables

t h 1 Integration
−
∫ k ⋅ dt = ∫ h
0 h0
2
dh

h
  1
k ⋅ t = 2 ⋅ h  2 Antiderivative
  h0
k2 2
h(t ) = ⋅ t + k ⋅ h0 ⋅ t + h0 Solution (parabola)
4

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-71-
 Determine the Cross Section Area

•Plot h(t) with assumed cross section AV

•Change AV until it fits

•Please note the change of water-jet contraction

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-72-
 Water Tank: MATLAB-Code
load [Link]
tm=FlowOutExperiment(:,1); %time [sec]
uS=FlowOutExperiment(:,2); %transducer voltage [V]
p=(uS-2)*1e4/8; %hydrostatic pressure [Pa]
rho=1000; %density of water[kg/m^3]
g=9.81; %acceleration of gravity [m/s^2]
hm=p/rho/g; %measured water level [m]
h0=0.308; %initial value [m]

d=8e-3; %assumed diameter of aperture of the ball valve [m]

AV=pi*d^2/4*0.8; %Cross sectional area of the ball valve [m^2]

AT=0.0028; %Cross sectional area of the tube [m^2]

k=-AV*sqrt(2*g)/AT; %abbreviation

tsim=0:0.01:13.5; %simulation time [sec]

hsim=k^2/4*tsim.^2+k*sqrt(h0)*tsim+h0; %simulated level [m]
plot(tm-0.75,hm,tsim,hsim,'r')
grid
xlabel('time [sec]')
ylabel('water level [m]')
title('Flow out experiment')
axis([-0.5 15 0 0.35])

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-73-
 Flow Out Experiment
0.35
measuremnt
simulation
0.3

0.25

0.2
level [m]

0.15

0.1

0.05

-0.05
-5 0 5 10 15 20
time [sec]

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-74-
 Heater Experiment
Halogen Temperature
lamp (20 W) Sensor LM 35

Banana jack
Voltage Supply

12 V

BNC socket ϑ 10 mV/°C

PWM
5V
0V

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-75-
 Circuit Diagram

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-76-
 Pulse Width Modulation

Problem: delivering power to a load (motor, heater, etc.) leads to power loss
and heat dissipation in the output stage.
A conventional linear output stage applies a continuous voltage to a load.

Solution: PWM
Just two system components are needed: a saw tooth wave generator and a
comparator

When should you use PWM? If the pulse period is much shorter than the
time-constant of the load, then PWM has a potential application. The time-
constants are defined, for example, by the thermal mass of a resistive heater
or the mechanical inertia of a motor. And it's the long time-constants that
average the pulses to a desired value set by the PWM's input voltage (low
pass character).

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-77-
 Pulse Width Modulation

uc
Control signal uc
uS
High frequency saw tooth uS

uc ≥ u S ⇒ u0 switched on
uc < u S ⇒ u0 switched off
t

u0
Average value of u0 is
proportional to control signal uc

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-78-
 Mean Voltage of PWM

Cycle duration
TP
u(t) ton
u0

_
u

t
T t on
1 P 1 ton
TP ∫0 ∫ u0 dt = u0 ⋅
u= u (t ) dt =
TP 0
TP

mean voltage is proportional to duty cycle ton

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-79-
 RMS of PWM

• Varying u continuously heat power is proportional to u2

• Root mean square urms

T t
1 P 2 1 on 2 t
u rms = ⋅ ∫ u dt = ⋅ ∫ u0 dt = on ⋅ u0
TP 0 TP 0 TP u2
u R P=
u 2
t on u ton 2 R
P= = rms
⋅ = ⋅ P0
0
R TP R TP

• With PWM: linear relationship between heating power and duty cycle ton !

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-80-
 Measuring the Step Response with NI 6024E-Board

•Step input from 40 % to 80 % PWM at t = 250 sec.

•PWM-frequency 25 Hz
•Sensor output voltage: 10 mV per °C, linear relationship
•sampling rate 1 kHz (sampling time 1 msec)
•input range ADC: -5 V to +5V
•ADC resolution: 12 bit ==> 10V/4096=2,5mV
40.15

•quantisation noise (or 0,25 °C) 40.1

40.05

40

temperature [°C]
39.95

39.9

39.85

39.8

39.75

39.7

39.65
44 46 48 50 52 54 56 58
time [sec]

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-81-
 Measured the Step Response

60

55

50

45
temperature [°C]

40

35

30

25

20
0 50 100 150 200 250 300 350 400 450 500
time [sec]

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-82-
 System Modeling
•No physical laws
•Here: assuming that the system behaves like a second order system due to the step
response ==> parametric system identification Θ K
G(s) = =
•Overdamped D>1 Y (T1 ⋅ s + 1) ⋅ (T2 ⋅ s + 1)
•Gain K K
=
•Time constants T1 and T2 T 2 ⋅ s2 + 2 ⋅ D ⋅T ⋅ s +1
•System input: duty cycle y; system output: temperature ϑ

Block Diagram
K 1
0 .4
T 1 .s+ 1 T 2 .s+ 1
Step input [-] Transfer Fcn Transfer Fcn 1 Scope

initial
42 .6 temperature
[°C ]

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-83-
 System Identification

real system

Step +
(…)2 Σ min
input -
model

K 1

T 1 .s+ 1 T 2 .s+ 1
Transfer Fcn Transfer Fcn 1

Strategy: minimising the sum of the squares of the errors

by varying the variables K , T1 , T2
MATLAB: fminsearch-function

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
[°C ]
Mechatronik u. Fahrzeugtechnologie
-84-
 Identification with MATLAB: Flow Chart

IdentifyHeater.m

call

SumSquaresErrors.m

call

[Link]

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-85-
 MATLAB Code
clear; close all
global t theta K T1 T2 k
load [Link]
theta=thetaStep40to80PWM25Hz(250000:500000); %rename [°C]
t=(0:length(theta)-1)*1e-3; %generate time vector [sec]
clear thetaStep40to80PWM25Hz
plot(t,theta)
grid
xlabel('time [sec]')
ylabel('temperature [°C]')
%initial values
K=15.1/0.4; %System gain [°C]
T1=2; %time constant [sec]
T2=1; %time constant [sec]
figure(2)
hold
k=0;
p=[T1 T2]
options = optimset('MaxFunEvals',60); %Maximum number of function evaluations allowed
p=fminsearch('SumSquaresErrors',p,options) %minimizes the sum of the squares of the
errors, parameter transfer
figure(2) Fakultät für Maschinenbau und Mechatronik

hold Modelling and Simulation of Dynamic Systems

Prof. Helmut Scherf -86-
Mechatronik u. Fahrzeugtechnologie
 MATLAB Code
function summe=SumSquaresErrors(p)
global t theta thetasim K T1 T2 k
T1=abs(p(1)) % time constant
T2=abs(p(2)) % time constant
if T1<0.01 %simulink crash preventing
T1=0.01
end
if T2<0.01
T2=0.01
end
sim('Step40percent') %start simulation
figure(1)
plot(t,theta,thetasim(:,1),thetasim(:,2),'r')
grid
xlabel('time [sec]')
ylabel('temperatur')
title('system identification')
%pause
q= theta-thetasim(:,2); %calculate errors
summe=sum((q.^2)) %calculate the sum of the squares of the errors
figure(2)
k=k+1; Fakultät für Maschinenbau und Mechatronik

plot(k,summe,'*') Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf

Mechatronik u. Fahrzeugtechnologie
-87-
 Results

system identification
60
p=
58
2 1
56

54
Temperature in °C

52
p=
50
24.3422 0.0127
48

46

44

42
0 50 100 150 200 250
time [sec]

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-88-
 Results
6
x 10
2

1.8

1.6

1.4
sum of squared errors

1.2

0.8

0.6

0.4

0.2

0
0 10 20 30 40 50 60 70
iteration steps [-]

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-89-
 Conclusions

•Identification tool works very well

•Assumption of first order system is sufficient
•Estimating system gain K is not necessary
•Always look upon a priori knowledge
•fminsearch is not a panacea

thuốc bách bệnh

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-90-
 Heating Experiment (Cooling Rib)

Aluminium rod
Electric
heating
device
LM 395 transistor
operating
in saturation mode

Temperature measurement points

Sensor: LM35 10 mV/°C

Lustre terminal
(power supply )

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-91-
 Set-up of the Differential Equation

dx L

PHeater Px Px + dx H

x Pconv B
We want to simulate the temperature T(x, t) which is a function of location x and time t.
P Heating power[W)]
Heater

Px , Px + dx Heating power due to conduction at location x, x+dx [W]

Pconv Heating power due to convection at location x [W]
B, H , L Width, hight, length [m]

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-92-
 Set-up of the Differential Equation

Thermal energy balance for the small mass element dm:

In words:
Incoming heat power minus outgoing heat power equals
the time rate of thermal energy change

∂Ethermal
Px − Px + dx − Pconv = (1)
∂t
Fourier’s law of heat conduction: ∂T
Heat conducted through cross sectional area Is Px = −k ⋅ Ac ⋅ (2)
proportional to the temperature gradient (minus ∂x
sign denotes heat flux always from high to low
temperature!) Radiation is neglected!
k Thermal conductivity [W/(mK)]
A c = B ⋅ H cross sectional area [m2]
∂T temperature gradient [K/m]
∂x
Fakultät für Maschinenbau und Mechatronik
Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-93-
 Set-up of the Differential Equation

∂P
Taylor series expansion: Px + dx = Px + ⋅ dx
∂x
∂ 2T
= Px − k ⋅ Ac ⋅ 2 ⋅ dx (3)
∂x
Thermal power loss due to convection is proportional to the heat
exchanging surface area and the temperature difference

Pconv = h ⋅ dA ⋅ (T − TAmb ) (4)

h Convection coefficient[W/(m2K)]
dA ≈ 2 ⋅ H ⋅ dx surface area [m2] (face surface neglected)
T Amb Ambient temperature [°C]

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-94-
 Set-up of the Differential Equation

time rate of thermal energy change in mass element dm

∂Ethermal ∂T
= c ⋅ dm ⋅
∂t ∂t
∂T
= c ⋅ ρ ⋅ dV ⋅
∂t
∂T
= c ⋅ ρ ⋅ B ⋅ H ⋅ dx ⋅ (5)
∂t
c Specific heat [J/(kg K)]
ρ
Density [kg/m3]
∂T
Derivative with respect to time [K/sec]
∂t

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-95-
 Set-up of the Differential Equation

Plug (5), (4) and (3) in (1):

∂ 2T ∂T
k ⋅ B ⋅ H ⋅ 2 ⋅ dx − h ⋅ 2 ⋅ H ⋅ dx ⋅ (T − TAmb ) = c ⋅ ρ ⋅ B ⋅ H ⋅ dx ⋅
∂x ∂t
Rearranging the equation:

∂T k ∂ 2T 2⋅h 2⋅h
= ⋅ 2 − ⋅T + ⋅ TAmb
∂t c ⋅ ρ ∂x c⋅ρ ⋅B c⋅ρ ⋅B
∂ 2T
= K1 ⋅ 2 − K 2 ⋅ T + K 2 ⋅ TAmb (6)
∂x
Second order partial differential equation

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-96-
 Initial Conditions, Boundary Conditions

Initial conditions:
At the beginning the whole aluminium rod has the same temperature TAmb

T ( x, t = 0) = TAmb
Boundary conditions:
Left boundary: we can calculate the temperature gradient at x=0 with Fourier’s law

∂T ∂T PHeater
PHeater = −k ⋅ B ⋅ H ⋅ ⇒ =− (7)
∂x x =0 ∂x x =0 k ⋅B⋅H
Right boundary: we have no heat transfer over right border at x=L

∂T
=0 (8)
∂x x =L

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-97-
 Simple Approach to Solving the Differential Equation
First step: Discretise the rod into 6 segments (length LS)

0 1 2 3 4 5

Second step: set-up the differential equation for each segment

Power balance segment 0 c ⋅ m0 ⋅ T&0

T0 − T1
PHeater 0 Pcond = k ⋅ Ac ⋅
LS

Pconv = h ⋅ A0 ⋅ (T0 − TAmb )

Fakultät für Maschinenbau und Mechatronik
Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-98-
 Power balance segment 1 c ⋅ m1 ⋅ T&1
T0 − T1 T1 − T2
Pcond = k ⋅ Ac ⋅ 1 Pcond = k ⋅ Ac ⋅
LS LS

Pconv = h ⋅ A1 ⋅ (T1 − TAmb )

Power balance segment 2
to 4 on the analogy of
segment 1
c ⋅ m5 ⋅ T&5
Pcond =
Power balance segment 5
T4 − T5 5
k ⋅ Ac ⋅
LS
Pconv = h ⋅ A5 ⋅ (T5 − TAmb )
Fakultät für Maschinenbau und Mechatronik
Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-99-
 Six Ordinary Differential Equations
T0 − T1
PHeater − k ⋅ Ac ⋅ − h ⋅ A0 ⋅ (T0 − TAmb ) = c ⋅ m0 ⋅ T&0
LS
T0 − T1 T1 − T2
k ⋅ Ac ⋅ − k ⋅ Ac ⋅ − h ⋅ A1 ⋅ (T1 − TAmb ) = c ⋅ m1 ⋅ T&1
LS LS
T1 − T2 T2 − T3
k ⋅ Ac ⋅ − k ⋅ Ac ⋅ − h ⋅ A2 ⋅ (T2 − TAmb ) = c ⋅ m2 ⋅ T&2
LS LS
T2 − T3 T3 − T4
k ⋅ Ac ⋅ − k ⋅ Ac ⋅ − h ⋅ A3 ⋅ (T3 − TAmb ) = c ⋅ m3 ⋅ T&3
LS LS
T3 − T4 T4 − T5
k ⋅ Ac ⋅ − k ⋅ Ac ⋅ − h ⋅ A4 ⋅ (T4 − TAmb ) = c ⋅ m4 ⋅ T&4
LS LS
T4 − T5
k ⋅ Ac ⋅ − h ⋅ A5 ⋅ (T5 − TAmb ) = c ⋅ m5 ⋅ T&5
LS
Fakultät für Maschinenbau und Mechatronik
Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-100-
 [T0]
PHeater PHeater [W ]
From 12
[W]
T0 [°C ] [T0]
[T1]
[T 1] T1 [°C ] Goto
From3
From1
segment 0
[T2]

From2
[T0] T0 [°C ]
[T3] Scope
From
T1 [°C ] [T 1]
From5
[T2] T2 [°C ] Goto 1
Simulink Block Diagram

[T4]
From 6
segment 1 From 11

[T5]
[T1] T1 [°C ]
From 15
From 4
T2 [°C ] [T 2]

[T3] T3 [°C ] Goto 2

From 8
segment 2

[T 2] T2 [°C ]

From7
T3 [°C ] [T3]

[T 4] T4 [°C ] Goto 3

From 10
segment 3

[T 3] T3 [°C ]

From 14 [T4] T4 [°C ] T5 [°C ] [T 5]

T4 [°C ] [T4]
From 9 Goto 4
[T 5] T5 [°C ] Goto 5

From 13 Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems
segment 4 [Link] 5
Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-101-
 Block Diagram: Interior of Segment 1

1 k*A/LS
T0 [°C]
Gain 3
T1' T1
1
2 k*A/LS 1/c/mS 1
s
T2 [°C] T1 [°C]
Gain 2 Gain Integrator

h*ASurf

Gain 1
Subtract
TAmb Constant

Initial condition: TAmb

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-102-
 Technical Data

Heating power PHeater = 7.35 W

4.8 V, 1.5 A (current limitation)
Thermal heat conductivity k = 174 W/(m ⋅ K)
Convection coefficient h = 13 W/(m 2 ⋅ K)
Specific heat capacity c = 1015 J/(kg ⋅ K)
Density ρ = 2700 kg/m 3
Width B = 2 mm
Height H = 20 mm
Length segment 0 LS 0 = 30 mm
Length segment 1 to 4 LS 1− 4 = 40 mm
Length segment 5 LS 5 = 45 mm
Ambient temperature TAmb = 26 °C

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-103-
 Simulation Result

120
T0
T1
110
T2
T3
100 T4
T5

90

80
temperature in °C

70

60

50

40

30

20
0 500 1000 1500
time in sec

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-104-
 Simulation Result

120

110

100

90

80
temperature in °C

70

60

50

40

30

20
0 50 100 150 200 250
distance in mm

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-105-
 Simulation and Measurement

100
T1
T2
90 T3
T4
T1sim
80 T2sim
T3sim
T4sim
70
temperature in °C

60

50

40

30

20
-200 0 200 400 600 800 1000 1200 1400 1600
time in sec

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-106-
 Numerical Solution of the Heat Differential Equation

Second order partial differential equation

∂T ∂ 2T
= K1 ⋅ 2 − K 2 ⋅ T + K 2 ⋅ TAmb
∂t ∂x

Initial conditions: T ( x, t = 0) = TAmb

∂T ∂T PHeater
Boundary condition: PHeater = −k ⋅ B ⋅ H ⋅ ⇒ =−
∂x x =0 ∂x x =0 k ⋅B⋅H

Finite difference approximation:

Replacing the partial differential quotient by difference quotient

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-107-
 Finite Difference Approximation
Forward difference approximation:

∂T T ( x, t + ∆t ) − T ( x, t ) Ti , j +1 − Ti , j
≈ = (9)
∂t ∆t ∆t
∂T T ( x + ∆x, t ) − T ( x, t )
≈
∂x ∆x
T ( x + ∆x, t ) − T ( x, t ) T ( x, t ) − T ( x − ∆x, t )
−
∂ 2T ∆x ∆x
≈ = ...
∂x 2
∆x
T ( x + ∆x, t ) − 2 ⋅ T ( x, t ) + T ( x − ∆x, t )
=
∆x 2
Ti +1, j − 2 ⋅ Ti , j + Ti −1, j
= (10)
∆x 2

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-108-
 Finite Difference Approximation
Plug (9) and (10) in (6):

Ti , j +1 − Ti , j Ti +1, j − 2 ⋅ Ti , j + Ti −1, j
= K1 ⋅ − K 2 ⋅ Ti , j + K 2 ⋅ TAmb
∆t ∆x 2

Solving this equation for Ti,j+1

K1 ⋅ ∆t
Ti , j +1 = ⋅ (Ti +1, j − 2 ⋅ Ti , j + Ti −1, j ) − K 2 ⋅ ∆t ⋅ (Ti , j − TAmb ) + Ti , j
∆x 2

Recursive algorithm with two loops

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-109-
 How the Algorithm works
1 i −1 i i +1 Nx x
1
∆x
j −1
j
j + 1X X
∆t
Nt XX
Calculation with
t boundary conditions

K1 ⋅ ∆t
Ti , j +1 = ⋅ (Ti +1, j − 2 ⋅ Ti , j + Ti −1, j ) − K 2 ⋅ ∆t ⋅ (Ti , j − TAmb ) + Ti , j
∆x 2

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-110-
 Considering the Boundary Conditions

∂T T2, j +1 − T1, j +1 PHeater

Left boundary: ≈ =−
∂x x =0 ∆x k ⋅B⋅H
PHeater ⋅ ∆x
Solving this equation for T1, j+1 T1, j +1 = + T2, j +1 X
k ⋅B⋅H

Right boundary: ∂T TN x , j +1 − TN x −1, j +1

≈ =0
∂x x=L ∆x
Solving this equation for TNx, jTN , j +1
x
= TN x −1, j +1 X
Fakultät für Maschinenbau und Mechatronik
Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-111-
 Simulation Result
Temperature Profile of heated Aluminium Rod
120

110

100

Steady state temperature

90

80
temperature in °C

70

60

50

40

30

20
0 50 100 150 200 250
Initial temperature distance in mm

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-112-
 Simulation Result

Temperature Profile of heated Aluminium Rod

120
1 cm
8 cm
110
16 cm
24 cm
100

90

80
temperature in °C

70

60

50

40

30

20
0 100 200 300 400 500 600 700 800
time in sec

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems 113
Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-113-
 Simulation Result

Temperature Profile of heated Aluminium Rod

0

60

100

50
200

40
300
time in sec

400 30

500 20

600
10

700

0 0.05 0.1 0.15 0.2

distance in mm

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-114-
 Steady State Solution of Heat Equation
∂T
Steady state means: rate of change equals zero: =0
∂t

∂T ∂ 2T
= K1 ⋅ 2 − K 2 ⋅ T + K 2 ⋅ TAmb = 0
∂t ∂x
∂ 2T
− K 3 ⋅ T = − K 3 ⋅ TAmb
∂x 2

Solution e.g. with MATLAB

syms K3 TAmb y x
dsolve('D2y-K3*y=-K3* TAmb ','x')
ans =

exp(K3^(1/2)*x)*C2+exp(-K3^(1/2)*x)*C1+TAmb

− K3 ⋅x K3 ⋅x
T ( x) = C1 ⋅ e + C2 ⋅ e + TAmb
Fakultät für Maschinenbau und Mechatronik
Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-115-
 Determination of arbitrary Constants
Boundary conditions:
∂T PHeater ∂T
=− =0
∂x x =0 k ⋅B⋅H ∂x x =L

∂T − K3 ⋅x K3 ⋅x
= −C1 ⋅ K 3 ⋅ e + C2 ⋅ K 3 ⋅ e
∂x

PHeater
x = 0 : −C1 ⋅ K 3 + C2 ⋅ K 3 = −
k ⋅B⋅H
− K 3 ⋅L K 3 ⋅L
x = L : −C1 ⋅ K 3 ⋅ e + C2 ⋅ K 3 ⋅ e =0
− K 3 ⋅L K 3 ⋅L K 3 ⋅2⋅ L
− C1 ⋅ e + C2 ⋅ e = 0 ⇒ C1 = C2 ⋅ e

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-116-
 Determination of arbitrary Constants

K 3 ⋅2⋅ L PHeater
− C2 ⋅ e ⋅ K 3 + C2 ⋅ K 3 = −
k ⋅B⋅H
− PHeater
⇒ C2 =
k ⋅ B ⋅ H ⋅ K3 ⋅ 1 − e ( K 3 ⋅2⋅ L
)
− K 3 ⋅L K 3 ⋅L
x = L : −C1 ⋅ K 3 ⋅ e + C2 ⋅ K 3 ⋅ e =0
− K 3 ⋅L K 3 ⋅L K 3 ⋅2⋅ L
− C1 ⋅ e + C2 ⋅ e = 0 ⇒ C1 = C2 ⋅ e

syms K3 TAmb y x k B H PHeater L

dsolve('D2y-K3*y=-K3* TAmb','Dy(0)=-PHeater/k/B/H','Dy(L)=0','x')

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-117-
 Plausibitlity Check
In the steady state the whole incoming heat power equals
the heat loss due to convection

L
PHeater = ∫ h ⋅ (T ( x) − TAmb ) ⋅ dA
0
L
= ∫ 2 ⋅ H ⋅ h ⋅ (T ( x) − TAmb ) ⋅ dx
0
Nx
≈ ∑ 2 ⋅ H ⋅ h ⋅ (Ti − TAmb ) ⋅ ∆x
i =1

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-118-
 Relay

contacts
Armature off
Armature on

coil

connections

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-119-
 Relay
n Number of windings [-]

i Current [A]
k spring
R Resistance [Ω]
l
x x0 L Inductance [H]

Air gap
H Magnetic field intensity [A/m]

Iron core B Magnetic flux density [T=Vs/m2]

n windings µ 0 = 4π ⋅10 −7 Vs/(Am)

permeability of free space,
absolute permeability

µr Permeability [-]

i u Φ Magnetic flux [Vs]

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-120-
 Electric Differential Equations

u = R ⋅ i + (L ( x ) ⋅ i )
Kirchhoff‘s voltage law:
d
dt
di dL
= R ⋅ i + L( x) ⋅ + ⋅i
dt dt
with

dL dL dx dL
= ⋅ = ⋅ x&
dt dx dt dx
Velocity of armature

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-121-
 Calculating the Magnetic Flux

Ampere‘s law: ⋅ n = ∫ H ⋅ ds ≈ H1 ⋅ l1 + H 2 ⋅ l2
i{ 1: iron
2: air
ampere-turn l Length of the field line

Φ
Neglecting leakage flux: B1 = B2 = B = = µ r ⋅ µ 0 ⋅ H1 = µ 0 ⋅ H 2
A
B1 B2
⇒ H1 = , H2 =
µ r ⋅ µ0 µ0

 l1 l2  Φ  l1 x 
i ⋅ n = B ⋅  +  = ⋅  + 
 µ r ⋅ µ0 µ0  A  µ r ⋅ µ0 µ0 
i ⋅ n ⋅ A ⋅ µr ⋅ µ0
Solving this equation for Φ: Φ =
l1 + µ r ⋅ x
Fakultät für Maschinenbau und Mechatronik
Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-122-
 Calculating the Inductance

Φ n 2 ⋅ A ⋅ µr ⋅ µ0
Inductance of a coil : L = n⋅ =
i l1 + µ r ⋅ x
Plug in Φ

Derive L with respect to x: dL n 2 ⋅ A ⋅ µ r2 ⋅ µ 0

=−
dx (l1 + µr ⋅ x )2

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-123-
 Calculating the Magnetic Force

1
Stored energy in a coil
W = ⋅ L ⋅ i2
2

Magnetic force F:

dW = F ⋅ dx
dW 1 dL 2 n 2 ⋅ A ⋅ µ r2 ⋅ µ0 2
⇒F= = ⋅ ⋅i = ⋅i
2 ⋅ (l1 + µ r ⋅ x )
2
dx 2 dx

Remark: i does not depend on x! See measured steady state values of current
With no air gap, with constant air gap, with variable air gap

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-124-
 Mechanical Differential Equation

− F ⋅ l + (Fpreload + k ⋅ ( x0 − x ))⋅ l = J ⋅ ϕ&&

Newton‘s law:
equilibrium of moments :

Kinematic relationship: x = l ⋅ ϕ ⇒ &x& = l ⋅ ϕ&&

− F ⋅ l + (Fpreload + k ⋅ ( x0 − x ))⋅ l = J ⋅
&x&
l
Fpreload + k ⋅ ( x0 − x )
l FSpring
x ϕ
Fpreload
F
0 x0 x
Fakultät für Maschinenbau und Mechatronik
Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-125-
 From 1 [F ]
Block Diagram x'
[xdot ] Goto 2
[F]

From 4
[I]
F_preload
x''
1 1 From 5
Constant l^2/J [x]
s s
u<0 [x]
Gain 1 Integrator
Integrator Goto
Product 2
Gain if FM <Fspring Product initial conditon :x0 From 7
k noacceleration
Fcn4 [xdot ]

u>0 From 8
x0 if gap is closed (x=0) Scope 2
== >xdot =0
Constant 1 [I]

From 6 Scope 5
From
[x]

u=12 V 12

Product 1
di /dt i
1
R
s
Integrator 2 Gain 2

Goto 1

[I] nW^2*A*myr ^2*my 0/2/(L1+myr*u(2))^2*u(1)^2 [F]

Fcn3 Goto 3

[x] -nW ^2*A *myr ^2*my0/(L1+myr *u(2))^2*u(3)*u(1)

From 2
dL /dx *xdot *i
[xdot ]
nW^2*A*myr *my 0/(L1+myr *u(2))
From 3
L(x)

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-126-
 Simulation Result

0.045

0.04

0.035
measurement
0.03 simulation

0.025
current in A

0.02

0.015

0.01

0.005

-0.005
-0.01 0 0.01 0.02 0.03 0.04 0.05
time in sec

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-127-
 MATLAB Code
%Technical Data from data sheet or measured or estimated
%---------------------------------------------------------
g=9.81; %constant of gravitation [m/s^2]
k=32e-3*g/1e-3; %(62-30) gr for 1 mm
F_preload=30e-3*g; %spring preload [N]
rho=7500; %density [kg/m^3]
l=12e-3; %length [m]
b=9e-3; %width [m]
d=1e-3; %thickness [m]
V=b*l*d; %Volume [m^3]
m=rho*V;
J=m/12*(l^2+d^2)+m*(l/2)^2; %inertia +parallel-axes theorem
[kg*m^2]
RShunt=22; %Shunt for current measure [Ohm]
R=288; %total resistance coil + shunt[Ohm]
my0=4e-7*pi; %absolute permeability [Vs/(Am)]
L1=70e-3; %length of field lines in iron [m]

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-128-
 MATLAB Code

x0=1e-3; %maximal air gap [m] and initial value integrator

aW=1e-3; %thickness of winding [m]
bW=18e-3; %length of coil [m]
dW=50e-6 %wire diameter [m]
nW=aW*bW/dW^2 %number of windings [-]

dC=7e-3; %core diameter [m]

A=pi*dC^2/4; %cross section of core [m^2]
myr=100; %Relative permeability [-]

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-129-
 MATLAB Code
load [Link]
load [Link]
load [Link]
figure
plot(withconstAirGap(:,1),withconstAirGap(:,2)/22,withnoAirGap(:,1),…
…withnoAirGap(:,2)/22,withvarAirGap(:,1),withvarAirGap(:,2)/22)
grid

pause
sim('Relay')
plot(withvarAirGap(:,1),withvarAirGap(:,2)/22,I(:,1),I(:,2),'r')
grid
xlabel('time in sec')
ylabel('current in A')
legend('measurement','simulation')

Fakultät für Maschinenbau und Mechatronik

Modelling and Simulation of Dynamic Systems Prof. Helmut Scherf
Mechatronik u. Fahrzeugtechnologie
-130-