Discrete time variational mechanics of multidomain systems: Applications to coupled electronic, hydraulic, and multibody systems

(1)

Discrete time variational mechanics of multidomain systems:

Applications to coupled electronic, hydraulic, and multibody systems

by

Tomas Sj¨ostr¨om

Master Thesis for the degree of

Master of Science in Engineering Physics

Supervisor: Claude Lacoursiere Examiner: Matrin Servin

November 23, 2012

Ume˚ a University Department of Physics

SE-901 87 UME˚ A

Sweden

(2)

(3)

Abstract

English

Today there exist few non-smooth multi-domain simulation tools using time- discretized Lagrangian mechanics for circuits.The main goal is to show that it is possible to use a semi-implicit, parameter free non-smooth variational time stepper to simulate the circuits with time-steps proportional to the system time scales.

This is demonstrated by implementing and performing extensive numerical tests for various types of electrical, mechanical and hydraulic components and demonstrate that the components are stable, with the correct behavior when the system is solved using a modified block pivot solver.

Simulation results shows that piecewise linear models are enough for the simple switching circuits in this thesis.

Svenska

Idag finns det f˚a simulatorer för icke-släta multidomän kretsar som bygger p˚a tidsdiskretisering av Lagranges ekvationer.

Huvudm˚alet är att visa att det är möjligt att använda en semi-implicit, parameter fri icke-slät diskret lösare för att simulera kretsar med tidssteg pro- portionella mot systemens tidsskalor.

Detta visas genom att implementera olika typer av elektriska, mekaniska och hydrauliska komponenter samt att visa att komponenterna är stabila och har rätt beteende när systemet simuleras av en modifierad block pivot lösare.

Simulerings resultaten visar att icke-släta Newton metoder med styckvis- linjära komponenter och komplementära villkor är tillräkligt för att simulera brytande komponponenter i de simulerande kretsarna.

(4)

Acknowledgments

First and for-most I want to thank my supervisor Claude Lacoursiere for helping me under a year will I worked on this Master Thesis. I also want to thank Algoryx Simulation AB for providing me with an interesting master thesis.

Last I want to thank all at UMIT Research Lab at Ume˚a University for letting me work there.

(5)

Contacts

Student

Name Tomas Sj¨ostr¨om

Adress Bobäcksvägen 12 92262 Tavelsjö Sweden Visit address MIT building, Campustorget 5, UMIT-lab Phone number 076 1470678

Email tomas sjostrom@hotmail.com Client

Name Algoryx Simulation AB

Adress Uminove Science Park, Box 7973SE – 907 19 Ume˚a Sweden Visit address Uminove Science Park, Box 7973SE – 907 19 Ume˚a Sweden Phone number 46-90-717090

Email contact@amgoryx.se Supervisor

Name Claude Lacoursiere

Adress Ume˚a University SE-901 87 Ume˚a Sweden Visit address MIT building, Campustorget 5, MA456 Phone number 46 706754242

Email claude@hpc2n.umu.se Examinator

Name Matrin Servin

Adress Ume˚a University SE-901 87 Ume˚a Sweden Visit address MIT building, Campustorget 5,MA456 Phone number 46 907866508

Email martin.servin@physics.umu.se

(6)

1. Introduction

Large vehicles consist of components from many different domains, a loader may have a mechanical scoop powered by hydraulic cylinders that is in turn controlled by an electric system. If the time-step used by the time stepper is longer then the time it takes to calculate the time-step, it’s possible to simulate the loader in real-time. It’s also important that the solution is stable with an error not diverting to infinity.

The main objective for this thesis is to evaluate an semi-implicit, parameter free non-smooth variational multi-domain time stepper that is build on work done for mechanical systems in Lacoursi`ere (2007). This is demonstrated by implementing and performing extensive numerical tests for various types of electrical, mechanical and hydraulic components and then demonstrate that the components are stable, with the correct behavior when the system is solved using a modified block pivot solver.

The thesis is divided into three parts. The first part describes the theory behind the project, the second part highlights special parts of the implementation and the last part contains results together withe proof of concept.

The theory part begins by explaining Lagrangian mechanics and constraints in chapter 3 with the equation of motion in equation 3.16. It counties in chapter 4 by showing that the different domains can be described using conservation of energy. Chapter 5 discretize the Lagrangian and shows the time stepper.

The chapter about linear complementarity problem shows how inequalities can be described by complementarity conditions. In chapter 7 the components are discretized and chapter 8 shows the dynamics for the different circuits.

The general principles behind the breadboard and modifications to the block pivot solver is in methods, chapter II.

The last part contains the results from the simulated circuits in chapter 10 and are then discussed in chapter 11. Chapter 12 together with chapter 13 contains the general conclusions from this master thesis.

The appendix’s contains the parameters from the simulations as well as an description of the principal sub matrices for an non-symmetric diode.

(12)

1.1. Problem Statement

The purpose of this thesis is to show, using extensive numerical testing that discrete time variational principle can be used to formulate and solve multi-domain systems. The goal is to systematic use differential algebraic equations(DAEs) with free variables to describe the system and avoid using co-simulation, but instead simulate the different domains in the same solver.

Solving DAE numerically often introduce numerical difficulties. The usual way is therefor to describe the circuits as systems of ordinary differential equations, ODEs. In this thesis this is avoided by using discrete-time variational mechanics.

The main reason for using discrete-time variational is to construct a solver that can simulate systems using time steps in the same order as the natural periods of oscillations of interest. This is unlike models like Spice requires unnecessary small time-step because of non-linearities represented by functions with sharp derivatives.

1.2. Related work

This is not the first time DAEs and the variational principle is used to described electrical circuits For instance Layton has a similar approach but unlike the description in this thesis, the equation of motion used by Laytons in (Layton, 1998, page 297) leads to singular matrices.

Ad hoc integration methods in Acary et al. (2011b) for non-smooth DAE’s includes free parameters, i.e., the θ method where they need to choose θ and don’t have a solid justification for it. Spice in Quarles et al. (2012) requires non-linear iterative solvers with small time-steps to resolve non-linear functions using smooth ODE’s. The Modelica solver in Modelica Association (2012) uses general high order, solver DAE solvers, while according to (Lacoursi`ere, 2007, chap 4) the solver in this thesis is a one step stable method tailored to for circuits. Marsden in Ober-Bl¨obaum et al. (2011) has Variational method for electrical circuits but not electronics. According to solving Modelica Association (2012) multi-domain systems with co-simulation using ODE integrators is gen- erally less stable then system DAE’s.

Many if not all these methods require clever preprocessing to eliminate the nodal matrix. For instance modified nodal analysis where the size of the solution matrix is reduced by making the smaller matrix more dense.

(13)

1.3. Original contributions

The theoretical work done in this thesis has been a collaboration between the author of this thesis Tomas Sj¨ostr¨om and the supervisor Claude Lacoursiere see table 1.1.

Tomas Sj¨ostr¨om Claude Lacoursiere

Work Section Work Section

Transformer 10.6 Resistor 7.2

None ideal motor 7.11.2 Inductor 7.5

Piston 10.12 Capacitance 7.6

The split circuit algorithm 9.2 Diode 7.7

Breadboard 9.1.1 Transistor 7.9

The expanded tikz library Whole report Op-amp 7.10

Motor 7.11.1

Pipe 7.13

Figure 1.1.: Shows the theoretic work done by Claude Lacoursiere

1.4. Limitations

Ideal components pose a problem as they can introduce undamped high fre- quencies. This means that if they operates at time scales shorter than the time-step won’t be damped. Two of the components, the op-amp and transistor introduce asymmetries in the complementarity problems. It is easier to solve symmetric positive define LCP problems. The proof to show that the asymmetries can be solved are not done. Only the time-step size has been compared between the different methods used in this thesis. The solution time has not been compared due to the solver used by this thesis is has not been optimized.

(14)

2. Definitions

Matrices in this thesis are denoted with capital letters iEX the impedance matrix is denoted by Z. Vector and scalars are denoted with small letters except for volume-rate that is denoted by Q, the distinction is apparent from the context.

Symbols

The following symbols are used in circuit diagrams Voltage Source Current Source

Resistor Inductor

Diode Motor

Transformer Opamp

NPN Transistor Mosfet Transistor

Ground Mass

Pump Check Valve

Spool Valve Piston

Relife Valve Reservoir

General circuit ^Z

(15)

Part I.

Theory

(16)

3. Lagrangian mechanics

3.1. Euler Lagrange equation

The theory in this master thesis is based on classical Lagrangian mechanics.

By using an energy based principle, the coupling between different domains, i.e. mechanical, electrical and hydraulic can be formulated in a unified way, and solved simultaneously. This is unlike methods i.e. Simulink that requires the use co-simulations using different time integrators and solvers for different domains.

Conservative physical systems can be described by the Lagrangian is given by

L^ps,s, t⁹ ^qT^ps,s, t⁹ ^qV^ps, t^q (3.1) where s are the generalized coordinates and s⁹ are the generalized velocities.

The action is defined as

S

»

L^ps,s, t⁹ ^qdt (3.2)

Equation 3.1 only accounts for conservative forces. In order to account for dissipative forces given by the Rayleigh dissipation function as defined in (Lacoursi`ere, 2007, page 68) to be

fℜ

Bℜ^ps,s⁹^q

Bs9 (3.3)

The energy dissipation rate for a Lagrangian without explicit depends given by (Lacoursi`ere, 2007, page 69) by

dE dt

Bℜ

Bs9s⁹ (3.4)

where E T V is the energy.

The derivation of the Euler Lagrange’s equation can be derived in different ways, using the Hamiltonian principle of least action or the energy principle as

(17)

in (Layton, 1998, page 290-297) and can be written as d

dt

BL

Bs9 BL

Bs

Bℜ^ps,s⁹^q

Bs9 (3.5)

3.2. Generalization

In Lacoursi`ere (2007) components are build on energy transference, storage and transformation which is common to mechanical, electrical and the hydraulic domain. Energy in a components is defined as

E

»

P^pt^qdt

»

pe^pt^qf^pt^{q q}dt (3.6) where E is the energy of the component, P is the power, e is the effort and f is the flow. For instance electrical circuits has the P IU. The systems are seen as ideal, energy don’t leave the system.

Due to the fact thats⁹and s is in general coordinates and can be interpreted as in table 3.1

Effort e Flow f Displacement s Momentum p

force F velocity v position x linear Momentum x

torque τ angular velocity ω angle Θ angular Momentum Θ

voltage U current I charge q flux linkage Λ

pressure P volume rate Q volume V Pressure momentum V temperature T entropy rate S entropy S ^pnone^q

Figure 3.1.: Shows correspondence between mechanical, electrical and hydraulic systems, from (Layton, 1998, page 276)

the energy postulate can be used in all domains.

3.2.1. Kinetic storage

The derivative of the momentum is defined as

p9e (3.7)

for example for a single particle with mass m the momentum is p mv and the effort is F mv.⁹

Replacing edt dp and using the same definition of kinetic storage as in (Layton, 1998, page 276) yields the kinetic energy part of the Lagrangian to

(18)

be

L^ps,s, t⁹ ^q

»

g^pf^qBf (3.8)

where g^pf^q is a function of the flow

Each component can then be derivatived with respect to the flow to yield the Euler Lagrange equation

d dt

BL^ps,s, t⁹ ^q

Bs9 g^pf^q (3.9)

3.2.2. Potential storage

The derivative of displacement flow as

s9f (3.10)

in mechanical systems vx⁹

Replacing f dt ds and using the same definition of kinetic storage as in (Layton, 1998, page 276) yields the potential energy part of the Largrantian to be

L^pq,q, t⁹ ^q

»

g^pf^qdf (3.11)

where g^pq^q is only depending on displacement.

Each component can then be derivatived respect to the displacement to yield the Euler Lagrange equation

BL^pq,q, t⁹ ^q

Bq g^pq^q (3.12)

3.2.3. Pure dissipators

Pure dissipator is by definition not conservative forces and can’t be described by the Lagrangian. A flow through a component exerts a dissipative effort on the system described in the same way as in (Layton, 1998, page 279) by a dissipative state function on the form

ℜ^pf^q

»

g^pf^qdf (3.13)

where ℜ is a Rayleigh dissipation function. For a resistor with resistance R and current I the state function becomes ℜ^pf^q ¹₂RI².

(19)

The derivative in the Euler Lagrange equation becomes ℜ^pf^q

Bf g^pf^q (3.14)

3.3. Constraints

Constrains are restrictions on the degrees freedoms in the system, in this thesis on the displacements, flows and on the effort. If a constraint is violated it applies a constraint force on the system that direct the system to a nonvolatile state, i.e. a ball penetrating a surface is force up to the surface by a non- penetration constraint.

In this thesis only holonomic(displacement dependent) g^ps, t^q ^¥0 and non- holonomic(displacement and flow dependent) a^ps, f, t^q^¥0 constraints are con- sidered. The f part of the Euler equation can in principle be seen as effort constraints with regularization in the form of mass. Pure effort constraints, without regularization are not allowed due to the fact that it is possible to put the constraints parallel to each other see figure 11.2.

The constraints can augmented to equation 3.1 as seen in (Lacoursi`ere, 2007, page 72-79) to

L

s,s, λ,⁹ λ, t⁹

L^ps,s, t⁹ ^q λg^ps,s⁹^q (3.15) subjected to Rayleigh function in equation 3.3.

The corresponding Euler Lagrange equation is derived in (Lacoursi`ere, 2007, page 72-79) to

d dt

BL^ps,s, t⁹ ^q

Bs9

BL^ps,s, t⁹ ^q

Bs G^TλA^Tα⁹ ^Bℜ^ps,s⁹^q

Bs9

d dt

BL

λ,λ, t⁹

B

λ9

BL

λ,λ, t⁹

Bλ g^ps,s, t⁹ ^q0

Bℜ^pα,α, t⁹ ^q

Bα9 a^ps,s, t⁹ ^q0

(3.16)

where the Jacobian matrix is defined as Gij

Bgi

Bqj and Aij

Bai

BIj. Note that the constraint force is given by Aα and Gλ. The magnitude of the force is given by λ ,α and it’s direction is given by the steepest descent (the Jacobian).

Equation 3.16 can be written

(20)

3.3.1. Charge-conservation

Coupling lumped components together is done using Kirchhoff’s current law or corresponding conservation of mass for hydraulics, i.e. that the current into a node must be the same as the current out from a node, if not the node will have a buildup of charge

The Kirchhoff’s current law can be written as

¸

i

a_jiq⁹_iI_s,j 0 (3.17)

where aji ^P^r1, 0, 1^sis the Jacobian matrix A also called nodal matrix and Is,j

is the sum of current source i at node j. The rows of the Jacobian corresponds to the nodes in the circuit and the coulombs corresponds to the currents. A component with one current connected to node one and two would have a Jacobian matrix on the form

A

1

(3.18) where the sign indicates if the current goes into the node or out from the node.

Equation 3.17 correspond to the holonomic conservation of charge constraint and is enforced by the current.

Equation 3.16 then becomes

d dt

BL^ps,s, t⁹ ^q

Bs9

BL^ps,s, t⁹ ^q

Bs A^TU ^Bℜ^ps,s⁹^q

Bs9

Aq⁹I_s

(3.19)

where f contains the current sources and the constraint force is the voltage U .

(21)

4. Multi-domain system

4.1. Mechanical systems

The flow variable for mechanical systems is as can be seen in table 3.1 velocity and the constraint effort is force. In this thesis only 1D kinetic stores in the form of 1D multi-body has been implemented to describe the mechanical world.

The kinetic energy becomes

T^pf^qT^pv^q 1

2mv² (4.1)

where m is the mass. The body can subjected to some potential V^ps^q , iEX gravity as well as some dissipative force g^pv^q. This yields the equation of motion3.5 to be on the form

mv⁹^BV

Bs g^pv^q (4.2)

and can then be augmented with other constraints, i.e. distance constraints as in equation 3.16

4.2. Electrical systems

The electrical systems are fundamental laws, Kirchhoff’s current(KCL) and Kirchhoff’s voltage law (KVL). KVL in electrical circuits is a result of conservation of energy and can be written as

¸

i

U_i 0 (4.3)

where Ui is the voltage drop over the component i. Equation states that 4.3 the sum of all the voltage sources in a loop must be equal to the sum of all the voltage drops over all loads in the same loop. In other case, energy will have been created or destroyed.

Coupling lumped components together is done using Kirchhoff’s current law.

The current into a node must be the same as the current out from a node, in other case the node will have a buildup of charge.

(22)

As seen in table 3.1 the effort for an electrical circuit is the voltage and the flow is the current I. In this thesis the electrical system consist of dissipators i.e.

resistors, kinetic storage’s i.e. inductors and potential storage’s i.e. capacitors.

The kinetic energy becomes

T^ps, f^qT^pI^q (4.4)

for inductors, the potential energy for capacitors becomes

V^ps^qV^pq^q (4.5)

and the dissipative state function becomes

ℜ^ps, f^qℜ^pq, I^q (4.6)

The equation of motion 3.5 becomes d

dt

BT^pI^q

BI

BV^pq^q

Bs

Bℜ^pq, I^q

BI (4.7)

which is correspondence to Kirchhoff’s voltage law due to that both laws ensure conservation of energy .

Equation 4.7 can then be augmented with constraints i.e. the conservation of charge as in equation 3.16.

4.3. Hydraulic systems

In correspondence to the electrical circuits in section 4.2 the hydraulically circuits are govern by the conservation of mass and the Bernoulli’s equation.

Bernoulli’s equation ensure conservation of energy is given in (Husain et al., 2009, page 31) as

P

ρ A²Q²

2 const (4.8)

where ρ is the density, A the area and the system has no high difference.

The connection between components are in the hydraulic case done with mass conservation. The mass of the flow into a node from one component must mach the mass of the flow out a node. In correspondence to equation 3.17 the conservation of mass becomes

¸

i

a_jiQ_iQ_s,j 0 (4.9)

(23)

where aji ^P^rA_i,0, Ai^sis the Jacobian matrix A, Qs,j is the sum of Q sources iat node j, A_i is the Cross section area of the components. A component with one flow connected to node 1 and 2 would have a Jacobian matrix on the form

A

(4.10) where the sign indicates if the fluid goes into the node or out from the node.

As seen in table 3.1 the effort for a hydraulic circuit is the pressure and the flow if the volume rate Q. In correspondence to the electrical circuits the hydraulic systems consist of dissipators by viscosity, kinetic storage’s through inertia in the fluid and potential storage’s in form of containers. The kinetic energy becomes

T^ps, f^qT^pQ^q (4.11)

for inductors, the potential energy for containers becomes

V^ps^qV^pV^q (4.12)

and the dissipative state function becomes

ℜ^ps, f^qℜ^pV, Q^q (4.13)

The equation of motion 3.5 becomes d

dt

BT^pQ^q

BI

BV^pV^q

BQ

Bℜ^pV, Q^q

BI (4.14)

which is correspondence to the Bernoulli’s equation due to that both laws ensure conservation of energy.

Equation 4.14 can then be augmented with constraints i.e. the conservation of mass from equation 4.9 as in equation 3.16.

(24)

5. Discrete time variational mechanics

5.1. Potential and kinetic storage

In this thesis the Lagrangian is discretized in the same way as in (Lacoursi`ere, 2007, page 48-49) to a sum of linear functions of short intervals ∆t

L^ps_k, s_{k 1},∆t^q

»h 0

L^ps,s, t⁹ ^qdt (5.1) where N is a constant, ss_{k 1} and s⁹ _∆t¹ ^ps_{k 1}s_k^q.

The discrete action is defined as

S^ps0, s1, .., s_N,∆t^q

N1

¸

k0

L^ps_k, s_{k 1},∆t^q∆t (5.2) The Hamiltonian principle of least action yields

BS^ps₀, s₁, .., s_N,∆t^q

Bs_k

BL^ps_k, s_{k 1},∆t^q

Bs_k

BL^ps_k1, s_k,∆t^q

Bs_k 0 (5.3) and the discrete time Lagrange equation

BL^ps_k, s_{k 1},∆t^q∆t

Bs^T_k

BL^ps_k1, s_k,∆t^q∆t

Bs^T_k 0 (5.4)

5.2. Pure dissipators

The derivative of the Rayleigh dissipation from equation 3.5 is done in the same way as in (Lacoursi`ere, 2007, page 66-68) to yield the discrete Euler Lagrange equation with dissipative force as

Bs^T_k

Bs^T_k ∆tfℜ^ps_k, f_{k 1}^q (5.5)

(25)

where fℜ is defined as in equation 3.3 of

Bs^T_k

∆tfℜ^ps_k, f_k^q (5.6) depending on which time-step the velocity is taken.

5.3. Component connection

The charge conservation from equation 3.17 is discretized at time-step k 1 and is enforce by the current. It is augmented to equation 5.5 as

BL^pq_k, q_{k 1},∆t^q∆t

Bq_k^T

BL^pq_k1, q_k,∆t^q∆t

Bq^T_k

∆tA^TU ∆tfℜ^ps_k, f_{k 1}^q ∆tU_e Af_{k 1}I_s

(5.7) where U is the constraint force force, Ueare external efforts and Isare external flows. Note that the constraint forces don’t depends on previous time-step.

5.4. Time stepping

For the case of inductors the kinetic energy has the form T^pt^q 1

2f^TLf (5.8)

where L is a block diagonal matrix with non-zeros one the diagonal for the kinetic components. For instance for a body with mass m T^pt^q ¹₂mv²

Potential storage can be can be written as V^pt^q 1

2s^TC¹s (5.9)

where C¹ is a block diagonal matrix with non-zeros on the diagonal for the potential components. For instances for a capacitor equation 5.9 would be V^pt^q _2C¹1q²

In this thesis only linearly dissipators are used this yields a Rayleigh dissipation function on the from

ℜ^pf^q 1

2f^TRf (5.10)

(26)

where R would contain the resistance on the diagonal in the same way as above.

Discrediting equation 5.8, 5.9 and 5.10 can be done in many different ways.

In this thesis they are chosen in order to yield a non-singular impedance matrix.

The kinetics storage are discretized with

f^pt^q s_{k 1}s_k

∆t (5.11)

The potential storage are discretized using the midpoint rule s^pt^q s_{k 1} s_k

2 (5.12)

The discrete Lagrangian becomes

L^ps_k, s_{k 1},∆t^q 1

4∆t²^ps_{k 1}s_k^q^TZ^ps_{k 1}s_k^q 1

8^ps_{k 1} s_k^q^TZ^ps_{k 1} s_k^q

(5.13)

The derivatived dissipation function becomes

fℜ^ps_k, f_{k 1}^qRf_{k 1} (5.14) Equation 5.7 divided by ∆t would then become

1

∆tL R ∆t

4 C¹

f_{k 1} A^TU Lf_kC¹

s_k∆t 4 f_k

e_e Af_{k 1} I_s

s_{k 1}s_k ∆tf_{k 1}

(5.15)

where e_e is eternal efforts and the solver first solves a matrix to get the flows and efforts then steps the displacements with an semi implicit Euler step.

Collecting the terms in front of f_{k 1} to an impedance matrices as

Z

1

∆tL 0 0

0 R 0

0 0 ^∆t₄ C¹

(5.16)

(27)

and defining the right hand side bLf_kC¹

s_k∆t 4 f_k

e_e (5.17)

as yields equation 5.15 to become

Z A^T

A 0

f_{k 1} e

b I_s

s_{k 1}s_k ∆tfk 1

(5.18)

where 0 is the zero matrix representation no relaxation of the constraint. The relaxation is not needed as the KCL is a linear constraint, the constraint is solved exactly and there is no need for relaxation.

Other constraints like the distance constraints may have relaxation with a diagonal matrix instead of zero matrix. This represents soften the constraint and is used in the original Spook solver to ensure that the solution matrix is a P -matrix. As can be seen in chapter 6, having no relaxation gives a Semi definite matrix requiring the algorithm from section 9.2 reduce the number of solutions to one when the circuit is split. The discretization is done in more depth for each component in chapter 7.

If you only remember one thing from this chapter remember the solution

matrix

Z A^T

A 0

(5.19)

(28)

6. The linear complementarity problem

The switching in this thesis is described using complementarity instead of fast changing function, i.e. the diode in section 7.7 has a complementarity condition on the form

I_s^¤I ^KU ^¥U_s (6.1)

which means that if I_s ^¤ I then U U_s and if U ^¤ U_s then I I_s. equation 6.1 is shown in figure 6.1

Us

I

Us UQ

IQ

Is

Figure 6.1.: Shows the ideal diodes IV curve

(29)

6.0.1. Inequality constraints

An Other example on complementarity are inequality constraints. The constraints have the form a^ps, f, t^q ^¥ 0, constraints that can be both active and inactive. A mechanical example would be a ball falling towards the floor with a non-penetrating constraint. When the ball falls towards the floor the non- penetration constraint is turn off. But when the ball penetrates the floor the constraint is activated and takes the ball to the surface.

In this thesis the inequality constraints will be used for describing diodes that blocks the current in one direction but not in the other. The constraint Euler Lagrangian is as described in (Lacoursi`ere, 2007, page 78-79) as

d dt

BL^ps,s, t⁹ ^q

Bs9

BL^ps,s, t⁹ ^q

Bs C^Tβ ^Bℜ^ps,s⁹^q

Bs9

0^¤c^ps,s, t⁹ ^q^Kβ^¥0

(6.2)

The^Kmeans that if the constraint force is negative the constraint is set inactive and if the constraint force is positive the constraint is active.

6.0.2. Impact stage

If the inequality constraint from section 6.0.1 was the only thing that stops the ball from falling through the floor in section 6.0.1 the ball would collide with the floor with non-zero incident velocity, and be forced to move away from the floor in a very short time interval. An impact model sets that time interval to zero.

This is imposed on the Lagrangian given in (Lacoursi`ere, 2007, page 167) as an LCP on the form

M G^T C_Imp^T

G 0 0

CImp 0 0

v λ ν

M v

Gv

w

0^¤ν^Kw^¥0

(6.3)

where M is the mass matrix, G is the Jacobian for other constraints, CImp is the Jacobian for the impact constraints, vis the velocity before the collision, v is the velocity after the collision, λ is the constraint force, w is a slack variable for the LCP. Notice that the n0 and the relaxations been removed.

When a collision is detected the the impact constraint activates and set the

(30)

velocity of the particle to zero. With the velocity set to zero the time-step before the impact the particle will never penetrate the wall.

6.1. Mixed LCP solver

The electrical and corresponding hydraulic circuits is described as an the linear complementarity problem(LCP) defined in the same way as in (Acary et al., 2011a, 51) as

wSxa

w^¯0, x^¯0, w^Tx0 (6.4)

where LCP is totally defined by the matrix S and vector a. The goal is to find a vector x that fulfills that Sxa^¯0 while x^¯0.

To ensure that equation 6.4 is solvable for any right hand side the S must according to (Lacoursi`ere, 2007, page 298) be a p-matrix, all principle sub matrices have non-zero determinant. Symmetric positive definite matrices are a subclass of p-matrices.

The solution matrix used in this thesis is positive semi definite p0-matrix with one or several solution, especially in the non-symmetric case with transistor and op amp the. To ensure that the LCP has one unique solution all principal sub matrix of S must have full rank. The LCP is solved with a block pivot algorithm. As can be seen in section 9.2 some of the principle sub matrices can have zero determinant. It’s possible to split the circuits into sub circuits with blocking diodes. This is taken care of in the same section by introduction artificial ground nodes.

A S on the form

S

Z A^T

A 0

(6.5) is solvable as long as the square nn matrix Z has full rank and the rectangular matrix mn A has full row rank. The complementary conditions used in this thesis is either on the currents corresponding to the Z matrix or on the voltages corresponding to the A. This means that either nodes or currents or both is set to inactive. In practices the equations are removed from the S and the corresponding variable is set to a fix value. If a row is removed the dimension Awill decrease but the rows will still be linearly independent and the LCP still solvable. On the other hand if a current is removed the system can split into

(31)

sub system. These is not the only restriction needed on the S to ensure that the LCP is solvable when the currents are remove.

All components in this thesis are build of Thevenin equivalent components consisting of bi-port devices with one current, i.e. the transistor has tree ports but is build by two diodes where each diode is connected to two ports. The most general form is one bi port device connected to two other circuits on the form

S

m 0 0 i j 0 k l 0 n

0 Z2 0 B1 B2 B3 0 0 0 0

0 0 Z3 0 0 0 B4 B5 B6 0

b A1 0 0 0 0 0 0 0 0

c A2 0 0 0 0 0 0 0 0

0 A3 0 0 0 0 0 0 0 0

d 0 A4 0 0 0 0 0 0 0

e 0 A5 0 0 0 0 0 0 0

0 0 A6 0 0 0 0 0 0 0

f 0 0 0 0 0 0 0 0 0

(6.6)

where b, c, d, e, f , i, k, l, m, n are real constants, m^¡0 ,Z are mass matrices, Aare a rectangular matrices with full row rank and Z B has full row rank. A is the nodal matrix and B is the corresponding matrix in KVL, not necessary symmetric.

When the bi-port device is connected to the circuit it can do it in six ways see figure 6.2

Appendix B in section 15 shows which constants that needs to be non-zero.

Note connection with entries not connected to any nodes are only allowed if the node is parallel to an other component also connected to these nodes. In other case it is possible to build a constraint loop which would require a new ground node. A loop would lead to a short circuit between the two nodes and in principle collapse them to one node. For instance the op-amp is working because the input side drives the out put side, but the output side don’t drive the input side. Note that it is possible to connected the input and output side but

In order for the principle sub matrices to have full rank the constants must fulfill table 6.3.

This shows that the components that are designed need to fulfill that the constants fulfill the requirement in table 6.3 as well as full rank in the S.

Note that this require the LCP to be solved by the modified block pivot non-smooth LCP sovler in order to handle the splitting of circuits.

(32)

Z3

Z2

1 2

3 4

5 6

Figure 6.2.: Shows the six different ways a component can be coupled between two circuits with impedance matrix Z2 and Z3

Connection type Required non-zero constants Required zero constants

1 f, n

2 b, d, i, k

3 b, f , i, n d, j

4 d, f , k, n e, l

5 b, c, i, j

6 d, e, k, l

Figure 6.3.: Shows the required values for the constant in equation 6.6 depending on connection type in figure 6.2. The mass is always m^¡0

(33)

7. Lumped components and their discretization

7.1. Creating component

The general method for creating electrical components are

* State the Kirchhoff’s volt law for the components. This is the equation that is to be solved for each components.

* Assume the KVL is the current derivative of the Lagrangian.

* Integrate the KVL with respect to current in order to get Lagrangian and Rayleigh function

* Augment corresponding constraints, i.e. charge conservation

* Discretization of the Lagrangian

* Derivative the Lagrangian with respect to charge at time k and voltage The hydraulic components subjected to lamina flow are done in the same way.

(34)

7.2. The resistor

The resistor can be seen as an ideal dissipator yielding a dissipative force on the system. This means that it can be described with a Rayleigh dissipation function.

Kirchhoff’s voltage law, KVL is given by

U RI (7.1)

The Energy dissipation is given by ℜ

»I 0

RIdI 1

2I^TRI (7.2)

A resistor has input ports, one in and one out node. This means that the Jacobian from section 3.3.1 becomes

A

1

(7.3) The resistors is discretized with I I_{k 1} as

ℜ R

2I_{k 1}² (7.4)

In order for it to have a term in the mass matrix.

The derivation of the Rayleigh dissipation function can be found in (Lacoursi`ere, 2007, page 68) and yields

Bℜ

Bq_k ∆tRI_{k 1} (7.5)

Bℜ

BU 0 (7.6)

Combining equation , , and 3.16 yields the discretized augmented Euler La- grange equation divide by ∆t to be

Zr 0

0 R

I_r I_{k 1}

A_r1 A_r2 1 1

U_i U_j

U_r 0

Ar1 1 A_r2 1

Ir

I_{k 1}

0 0

(7.7)

where Zr, Ari, Is, Urbelongs to the rest of the circuit and Ui, Uj is the voltage

(35)

at the nodes connecting the resistor to the circuit.

(36)

7.3. The ideal voltage source

An ideal voltage source is an infinite external source of energy imposing a voltage difference UExtbetween in between two nodes. In reality this is not the case, but there is also a small resistance R in a battery reducing the voltage to U₀. This gives a KVL as

RIU0 UExt (7.8)

The battery is connected to the circuit in the same way as a resistor with a Jacobian matrix as in equation . This yields the calculations for the battery to be the same as in section 7.2 with the different that the resistor now has an imposed external voltage on the rightsized. In correspondence to equation 7.9 the discretized augmented Euler Lagrange equation divide by ∆t for the battery becomes

Zr 0

0 R

I_r I_{k 1}

A_r1 A_r2 1 1

U_i U_j

U_r U_Ext

Ar1 1 A_r2 1

Ir

I_{k 1}

0 0

(7.9)

where Zr, Ari,Is ,Ur belongs to the rest of the circuit and Ui, Uj is the voltage at the nodes connecting the resistor to the circuit.

(37)

7.4. The ideal current source

A current source is an infinite supply of external voltage to two existing nodes in the circuit. In order to ensure that the current source is added to existing nodes it is attached in parallel to a resistor with resistance R.

This yields the same KVL as for the resistor but a different in the KCL. The current from the current source IExt will flow into one of the resistors nodes and out from the other the KCL form equation 3.17 becomes

¸

a_jiq⁹I_j 0

¸

a_kiq⁹ I_k0 (7.10)

As can be seen in equation 5.7 this corresponds to adding IExt on the right hand side of the current constraint. In correspondence to equation 7.9 the discretized augmented Euler Lagrange equation divide by ∆t for the battery becomes

Zr 0

0 R

I_r I_{k 1}

A_r1 A_r2 1 1

U_i U_j

U_r 0

Ar1 1 A_r2 1

Ir

I_{k 1}

IExt

(7.11)

where Zr, Ari,Is ,Ur belongs to the rest of the circuit and Ui, Uj is the voltage at the nodes connecting the resistor to the circuit.

(38)

7.5. The inductor

Inductors can be seen as kinetic storages. The voltage over the inductor is given by

U LdI

dt Ld²q

dt² (7.12)

The voltage U can be seen as the force created by the second derivative of the charge. The associated flux linkage, momentum is given as

Λ

»T 0

LdI

dtdtLI (7.13)

where t is an arbitrary time. The kinetic energy is defined as in (Layton, 1998, page 276-278) to be

T

»I 0

LIdI 1

2I^TLI (7.14)

The inductor has the same type of connections as a resistor, one in and one out node. This means that the jacobian from section 3.3.1 becomes

A

1

(7.15) The current is I ^dq_dt and is discretized to I ^q^k_∆t¹^q^k, as the equation is not stiff and will have a mass. The discrete Lagrangian at time k becomes

L^pq_k, q_k,∆t^q∆t 1

2∆t^pq_{k 1}q_k^q^TL^pq_{k 1}q_k^q (7.16) The discrete augmented Euler Lagrangian, divided by ∆t becomes

Zr 0 0 _∆t¹ L

I_r I_{k 1}

A_r1 A_r2 1 1

U_i U_j

U_r

1

∆tLI_k

A_r1 1 A_r2 1

I_r I_{k 1}

0 0

(7.17)

where Zr, A_ri, U_r,I_s belongs to the rest of the circuit and U_i, U_j is the voltage at the nodes connecting the inductor to the circuit.

(39)

7.6. The capacitor

The Capacitors can be seen as ideal potential storage with a voltage-drop given by

U q

C (7.18)

According to is the potential energy is given by (Layton, 1998, page 277-279) V

»q 0

Udq

»q 0

q

Cdq q²

2C (7.19)

The capacitor is a two device with an in and an out port. This gives a Jacobian to be

A

1 1

(7.20) In-order for the capacitor to have a term in the mass matrix the charge is desecrated with the midpoint rule as q^q^k ¹₂ ^q^k

The discrete Lagrangian becomes L^pq_k, q_k,∆t^q∆t h

8^pq_{k 1} q_k^q^T C¹^pq_{k 1} q_k^q (7.21) Using that Ik^pq_{k 1}q_k^qthe discrete augmented Euler Lagrangian, divided by ∆t can be written as

Z_r 0

0 ^∆t₄ C¹

Ir

I_{k 1}

Ar1 Ar2

1 1

Ui

U_j

U_r

1

C q_k ^∆t₄ I_k

A_r1 1 Ar2 1

I_r I_{k 1}

0 0

(7.22)

where Zr, Ari, Ur,Isbelongs to the rest of the circuit and UiU_j is the voltage drop over the capacitor.

Discrete time variational mechanics of multidomain systems: Applications to coupled electronic, hydraulic, and multibody systems