Modelling and Control of Fast-Switching Solenoid Direct Injection Valves Using a New Magnetics Library

This paper deals with the special challenges in modelling and controlling fast-switching solenoid valves with Modelica. For this, a solenoid injector for application in direct injection combustion engines with switching times around 100 μs is being used as an example. The occurring eddy current losses as well as the local saturation phenomena based on magnetic field displacement require for detailed network models of the magnetic domain. Therefore, a new Modelica magnetics library based on a consistent systems modelling analogy is being proposed and implemented, which increases solvability of the injector models compared to the Modelica Standard Library. Additionally, different one-dimensional contact mechanic models for representation of the bouncing behaviour of switching-type solenoids are evaluated. The complete electro-magnetomechanical model of the injector is then used for synthesis of novel closed-loop control schemes found by model inversion and parameter optimisation.


Introduction
The past decades have been shaped by development of many alternative vehicle powering concepts, however the combustion engine still has its benefits when it comes to, e.g., refuelling or power density of the storage medium. On the other hand combustion engines are still an active field of research due to possible improvements in terms of pollution and carbon dioxide emissions.
Most of today's engines are based on the direct injection principle, by which the combustion can be controlled more precise due to possible pre-and post-injections as well as higher variability with the main injections per rotation (Bauer et al., 2004, pp. 146 -147). The demands on injectors are thereby increasing; in particular the minimum switching times have to lie below 1 ms, while the injection pressures the valve has to deal with are also growing (Ströhla, 2012, p. 13). Due to higher manufacturing costs the use of piezo-actuated valves instead of solenoids is not an option for low-priced cars or engines.
These demands on solenoid injectors are met with a two-stage movement of the armature, which allows the latter to accelerate without opening the valve in the first stage. Once the armature has traversed this free stroke space and has accumulated some momentum, the valve seat is then opened by a needle, which gets picked up by the already moving armature (Denk, 2018, p. 17). Figure 1 shows such an injector with two stages of movement for the gasoline direct injection, on which this paper is based. Figure 2 illustrates the two-stage movement.
In order to minimise switching times and to overcome the time delaying influences of the coil inductance as well as the eddy currents, an electric excitation with a boost phase is required, as depicted in Figure 3. The boosting voltage of 65 V is used to "pump" huge amounts of energy within a short period into the magnetic system (Bauer et al., 2004, p. 153). After this, there are two phases during which the electric current is held at specified levels to provide different magnetic reluctance forces for catching and holding of the armature, respectively (Denk, 2018, p. 27).
The relatively large time delaying influence of the electromagnetic domain, compared to the switching times, limits the controllability of armature and needle positions and therefore of the injection process itself. The applied excitation scheme also leads to an overshoot in terms of armature velocity, which results in strong bouncing at the stop plate. This not only increases noise, but also reduces life span of the components due to mechanical wear. In addition, wear is also a major source for changing injection amounts during life time (Rauer et al., 2019). Another demand for the injection of smallest amounts of gasoline is the operation in so-called ballistic mode, in which armature and needle do not reach the stop plate (Denk, 2018, pp. 18 -19). While this allows the engine to run in alternative combustion modes, the sensitivity of the output spray based on deviations in manufacturing and assembling or wear is still a problem.
Therefore, this paper addresses the challenges when modelling the electric, magnetic and mechanic domains of an injector in order to use this model for synthesis of novel control schemes. Fluid mechanics is neglected here, because lumped-parameter models are not suited for the complex behaviour of the injection fluid inside the valve, where shock waves and combustion pressures have to be taken into account. This domain could be included later by connecting CFD software via FMI, e.g.
2 Modelling the magnetic domain 2.1 Magnetic network analogies For system simulation purposes it is common to transfer the usually as partial differential equations formulated physics into a network of lumped elements. As connection variables there can be chosen multiple sets of pairs consisting of a potential and a flow variable. For the magnetic domain usually the magnetic flux Φ gets utilised as the flow variable, for which reason the magnetic reluctances R m are understood as "magnetic resistors". This is due to the fact, that component equations connecting potential and flow variable derivative-free with each other are commonly assumed to show dissipative behaviour like e.g. the electric resistor (Zimmer and Cellier, 2006). However, this does not hold for magnetic reluctances, as they merely store energy in a capacitive way. To address this, the magnetic flow variable has to be changed to the time derivative of the magnetic fluxΦ, which leads to the reluctance component equatioṅ with the magnetic potential V m and the reluctance capacity C m = 1/R m (Littmann and Schiedeck, 2008). In demarcation against the classic analogy we want to call this Table 1. Analogies describing magnetism by lumped elements.

Classic analogy Consistent analogy
principle the consistent analogy throughout this paper, as the (time-derived) potential and flow variables in the component equations indicate the correct energy-related behaviour. Additionally, the product of potential and flow variable is always a power value. Similarly to the reluctance, we can also rewrite the equations of eddy current elements -expressed in the magnetic rather than the electric domain -by introducing the parameter L m ≡ R * m ≡ 1/R ec as summarised in Table 1, with R ec being the electric resistance of the eddy current path. A physically inductive element in the magnetic domain does not exist.
As a side effect, the coupling equations between electric and magnetic domain are also simplified, since both AMPÈRE's law (Kallenbach et al., 2017, p. 18) and FARADAY's law (Kallenbach et al., 2017, p. 30) are now defined via the pure flow variable of the counterpart domain (I andΦ) multiplied with the winding number N, respectively. Here, the electric domain is represented by potential U and flow I. This coupling principle can be observed in energy transformations of other physical domains and helps in teaching and understanding generalised mechatronic systems modelling (Grabow, 2013;Littmann and Schiedeck, 2008). The classic analogy is widely used for network-based simulation purposes as well as for engineering-level rough design tasks (Birli et al., 2003;Ströhla, 2012). Also the Modelica.Magnetics.FluxTubes library bases upon this type of analogy (Bödrich and Roschke, 2005;Bödrich, 2008), which was extended by models for covering magnetic hysteresis phenomena (Ziske and Bödrich, 2012). FRANKE et al. tested the library with nonlinear magnetic networks, but did not include eddy current models (Franke, 2012). However, the FluxTubes library has not been tested yet on the literature with more complex magnetic networks consisting of eddy current and magnetic field displacement models, which are inevitable for all types of fast switching solenoids.
The present work has tried to apply the MSL FluxTubes library to medium sized magnetic networks with a total of 50 to 60 elements, containing nonlinear ferromagnetic material behaviour as well as eddy current elements without the use of the PREISACH or TELLINEN hysteresis models. Dymola and OpenModelica failed either in determining the causality or in solving the resulting system of equations, depending on software version and magnetic network structure, until the eddy current elements were excluded from the network 1 . A reason for that could not finally be determined despite the fact, that the use of eddy current elements seems to be one trigger. Therefore, in the following section an approach based on the consistent magnetic analogy is being developed, with the aim of higher solvability when applied to complex eddy current enhanced networks.

A new Modelica magnetics library based
on the consistent magnetic analogy Implementation The connector class of the new, FluxTubes_PhiD-called library basically has to use the time derivative of the magnetic fluxΦ as the new flow variable. Even with this change of connection variables, the KIRCHHOFF current law (all flows sum to zero at every node) still applies, what can be shown by differentiation: On the other hand, the changed formulation of the reluctances introduces one new state for every reluctance element, which now in general needs an additional initial value for integration and hence complicates initialisation. When considering non-excited magnetic circuits with Θ = 0 in steady state, it becomes clear that all magnetic potential differences ∆V m,j must be initialised with zero, as long as the magnetic circuit is unpolarised (Kallenbach et al., 2017, pp. 61 -63). In every other nontrivial case the initialisation problem is not well defined exclusively via the consistent analogy, as the conditions for initialisation get lost due to the derivation of the magnetic flux in the connectors. An additional use of the classic analogy only for initialisation could be promising, however, this is left to future research. Therefore, simulations in this paper all start with a non-excited magnetic circuit.
Similar to this, the linked magnetic flux Ψ can not be calculated out of the momentarily total flux Φ by but instead needs to be integrated from the flux derivative: Again, in steady state initialisation the initial value for the linked flux is assumed to be zero. Since the linked flux does only serve as an output variable and is not being used in other equations, this integration over the simulation period is not causing precision losses elsewhere. The same holds for the flux density B of each reluctance pathway l, which needs to be calculated via the magnetic permeability µ = µ 0 µ r out of the field strength H: As the consistent analogy does not swap potential and flow variable with each other, the topology of the network remains the same (Littmann and Schiedeck, 2008). Therefore in a third step only the component equations have to be reformulated according to Table 1. It should be noted though that in case of ferromagnetic materials, starting from HOPKINSON's law, the time derivative of the magnetic permeance G m ≡ 1/R m can not be neglected in general: However, in case of linear, structure-invariant reluctances of constant pathway length the termĠ m V m becomes zero and hence the necessary iteration depth for the numerical solution decreases. Here, the Evaluate = true flag should be provided for efficient code generation, since the assumptionĠ m V m = 0 does only hold for special parameterisation cases.
If, in addition, the magnetic permeability correlation µ r (H) for ferromagnetic reluctances is inserted into the term ∂ G m /∂ Φ and solved manually instead of being implemented into a causalised function, Equation (8) can be converted into the less implicit shapė where C m is nonlinear in V m , but the time derivativė G m is disappearing. Herewith, the numerical integration of ferromagnetic reluctances can be accelerated due to fewer numerical root searches, compared to Equation (8). This has been implemented as a separate library called FluxTubes_PhiD_fast with the drawback of a noninterchangeable permeability function.
Validation The Modelica libraries described above have been tested against the MSL FluxTubes implementation as well as a commercial simulation program called SESAM (Birli et al., 2003) in order to determine possible errors in the new implementation. This has been done by comparing the calculated reluctance forces for the example network SimpleSolenoid of the FluxTubes library (Bödrich, 2008). Apart from numerical aberrations, no systematic deviations in the network calculations have been noticed. In the simulations described later, also no drift due to using the flux derivative has been noticed.  In a second step, several different implementations of magnetic libraries underwent an analysis w.r.t. the solvability of complex magnetic networks. For this, two different circuits -one with consideration of eddy currents and the other one without -were prepared based on the lumped-parameter model of the injector presented in Section 2.3. These networks were tested with quasistationary boundary conditions (fixed armature movement) and in dynamic configuration (free moving armature). Additionally to the consistent-and classic-based libraries described above, other implementation types of the same analogies were tested, which differ only in the way the magnetic permeability gets computed. This is due to the fact, that the iterative calculation of the nonlinear material characteristic µ r (H) or µ r (B) has been observed to be another crucial point for solving the network. Table 2 summarises all herewith tested libraries with a short explanation of the implementation details and provides the calculation timings, if solving the network was possible 2 .
First of all, the two implementations of the consistent analogy described above -FluxTubes_PhiD and FluxTubes_PhiD_fast -lead to a solvable system of equations in all test cases and with both tested programs. The variant based on Equation (9) performs moderately better with speed advantages between 5 % and 30 %.
The MSL libraries based on the classic analogy seem to be harder to solve than the consistent analogy; the simulation process fails in most cases, either due to problems with causalisation or simulation. The few test cases, in which the simulation succeeds without complications, allow the conclusion that the consistent analogy is considerably more time-intensive to solve, which can be explained by the additional state variablesV m being introduced in the reluctances. Separate comparative simulations of small networks indicate the consistent analogy implementation performing almost an order of magnitude slower.
Inside the consistent analogy based libraries the permeability µ r is expressed via the magnetic field strength H instead of the flux density B. Additionally, the correlation µ r (H) is provided as an equation rather than a causalised function. In order to check, whether the definition of the permeance has an influence on the solvability, the two additional libraries FluxTubes_MuRAsEq and FluxTubes_MuROfH based on the classic analogy have been tested as well (c.f. Table 2). Herewith the solvability problems are still present, so that this approach seems also not promising. Furthermore, a separate test with the permeance expressed by an interpolation table resulted even worse, which is understandable due to more complicated differentiation of the look-up table output.
Hence, for the tested compiler versions of Dymola and OpenModelica the implementations based on the consistent analogy are assumed to be stable enough for modelling of complex magnetic circuits. The following simulations are all based on these libraries. It should be noted though, that Table 2 delivers only a snap-shot in terms of tested software, compiler versions and the specific magnetic network. E.g., in OpenModelica 1.14.2 occasionally some problems with the FluxTubes_PhiD_fast library are occurring, whereas the library FluxTubes_PhiD shows no problems.

Deriving the lumped-parameter network
Field displacement In ferromagnetic materials usually eddy currents are evoked proportionally to the magnetic flux derivativeΦ. Moreover, in axially permeated solids of revolution a field displacement effect similar to the skin effect can be observed, which leads to a magnetic field "invading" from the outermost to the innermost radius. This is caused by induced rotating electric currents, which again excite a magnetic potential ∆V m against the initial direction. Herewith, the radially inhomogeneous field with local saturations of the ferromagnetic material leads to  varying reluctance forces of the armature.
For this effect STRÖHLA found a substitute network with cascading reluctance and eddy current elements, as depicted in Figure 4 (Ströhla, 2002). Each pair of reluctance R m and "magnetic inductance" L m represents one shell of the body with constant cross sectional area (Kallenbach et al., 2017, p. 170). In order to prevent discontinuous coil currents, which can not be observed in reality, a small additional reluctance R m,0 should be inserted in a way, that the overall reluctance magnitude stays the same. In reality, this fact is given by the variety of stray fields over non-conductive materials of the actuator. Depending on the wall thickness of the hollow cylinder a discretisation of n ∈ [3 , 5] is sufficient for most applications, the result of which is demonstrated in Figure 5. The described model is implemented in the new libraries with variable discretisation depth n and propagation of the model parameters to the inner components, so that only the electric resistivity has to be provided additionally to regular reluctance elements.
Magnetic network With all basic elements and libraries being implemented, the lumped-parameter network can be derived from the injector geometry, the result of which is depicted in Figure 6. Like with all pot-shaped solenoids the reluctance pathways can be reduced to a two-dimensional, rotationally symmetric form, where the use of radially and axially permeated reluctances is sufficient. In case of highly saturated solenoids, stray fields over the coil area have to be taken into account (c.f. white elements in Figure 6), since the saturated ferromagnetic permeability converges to the vacuum permeability µ 0 . The inclusion of different stray pathways usually needs an iterative procedure, while for the coil unified approaches are known (Ströhla, 2012, pp. 35 -42). The more saturated parts of the solenoid exist, the less accurate the lumped-parameter approach will be.
The radially permeated reluctance elements in Figure 6 are modelled via an additional eddy current element, whose current pathway is assumed to be along the circumference of the mid radius. The majority of the remaining ferromagnetic pathways are represented by the field displacement models. However, one major drawback of this network topology is the fact, that the field displacement models of outer hollow cylinders do not influence the inner ones. Therefore, the total eddy current losses as well as the local saturation effects are presumably estimated too small in the model.
Validation In order to evaluate the accuracy of the magnetic network yet without other physical domains, a stationary comparison between a finite element model and the Modelica model is shown in Figure 7. Starting with the linked magnetic flux Ψ , the influence of saturations in the injector is obvious, as the model assumes the magnetic excitation as too small. This is a common aberration when transforming partial differential equations into a system of ordinary differential equations accompanied with lumped parameters in a network 3 , since the simplification of magnetic pathways results in underestimating the total permeance of the system, as long as no pathways are mapped redundantly. On the other hand, the progression of the linked magnetic flux Ψ w.r.t. the air gap δ m is well described by the Modelica model.
The comparison of the reluctance force F m shows a similar result, but with higher divergence between the depicted electric excitations. Here, the discretisation of the working air gaps (c.f. Figure 6) shows a high sensitivity w.r.t. the shape of the reluctance force curve. The serial combination of each air gap model with a short ferromagnetic pathway for modelling local saturation effects is also fundamental.
Further increased accuracy could be received by optimising the network manually based on the magnetostatic results or by conducting a parameter identification with the help of numerical optimisation algorithms (Mühlenhoff et al., 2016). However, this would be beyond the scope of the network modelling approach, since a finite element model is more suitable in cases, where calculation times aren't a problem. Since the overall model aberrations in Figure 7 lie within 5 % to 10 %, the network based approach is sufficient for the desired purpose of optimising the injector control. search. The movements of armature and valve-needle can easily be described by one-dimensional and translational equations of motion and do not require further investigation. On the other hand, contact problems in mechanics are known for their numerical stiffness and applicationspecific implementations (Tiller, 2001, pp. 162 -166). Therefore, in this chapter different discrete and continuous contact models will be compared in terms of modelling effort and integration performance. The injector to be modelled has a total amount of four end stops, with two of them acting in-between the bodies. Therefore, both masses are restricted in terms of position. Figure 8 shows the topology of the contacts neglecting other mechanical elements, such as acting forces.

Restitution-based contact models
One approach for modelling kinetic contacts is to reinitialise the state variables for the velocity v, whenever a previously defined boolean contact condition gets fulfilled. The principle of conservation of linear momentum is then used in order to determine the velocities v after contact. Therefore, this approach is limited to simultaneous contact of only two masses. If a partial loss of energy is assumed by introducing a restitution coefficient c r ∈ [0 , 1], the reinitialisation can be described by and vice versa (Gross et al., 2004, p. 93).  This hybrid model is known for its ZENO-behaviour, which leads to infinite high computation effort to approximate a stationary contact (Heymann et al., 2005), since the static contact constraint of resting masses is not explicitly considered. Another numerical problem arises when bodies begin to fall into each other, caused by the reinitialised velocities being too small. This is often treated by detecting invalid velocities after reinitialisation and assigning the corresponding acceleration to zero for simulating "sticking" behaviour (Tiller, 2001, p. 164).
For the injector implementation, one has to use a state graph with a total of 19 each-exclusive boolean states and 41 transition conditions, in order to cover all possible combinations of contact situations. The state machine has been implemented by hand into the mechanics model without the use of specific state graph libraries. Though, the resulting model performed poor w.r.t. integration performance, as Table 4 shows 4 . In order to cover all these drawbacks and to enable reuse of the models in an objectoriented modelling fashion, the following section investigates different non-phenomenological models.

Force-based contact models
KELVIN-VOIGT models Force-based contact models share the basic principle of a contact force F c being composed out of a capacitive and a dissipative part as functions of the indentation δ . The capacitive element can be modelled by a stiff linear or nonlinear spring force F s , which represents the deformation of both contact surfaces. The dissipative part F d is needed in order to consider energy losses. In case of a constant damping coefficient d and stiffness c, we get the linear KELVIN-VOIGT contact model (Machado et al., 2012): Because contact forces are intended to be always negative or zero -as long as no adhesion of the colliding bodies should occur -we can modify Equation (11) to which leads to substantially higher contact exit velocities. Figure 9 shows these resulting differences in the force-4 DASSL was used as integration method with a tolerance of 10 −5 . The failed simulations exhibited convergence problems at contact time points. The given simulation timings were averaged out of 3 trials. indentation-diagram. Nevertheless, both variants still exhibit the problem of discontinuous forces at contact start, resulting from the dissipative term in Equation (11), which is assumed to be unphysical (Machado et al., 2012). The ElastoGap model of the MSL handles this by multiplying the damping term dδ with the indentation δ (c.f. Figure 9).
However, the primary problem of the KELVIN-VOIGTbased models is the parameterisation of the damping coefficient d. While the stiffness rate c can be guessed reasonably via the surface geometries, there are no accurate approaches known for estimation of the damping constant as a function of the restitution coefficient c r , which can be measured more easily and is less dependent on the initial impact velocity. In contrast, the damping constant d strongly depends on the relative velocities and therefore the model can not be deployed in varying boundary conditions with reasonable validity.
Hysteresis damping models This parameter identification problem can be handled by relating the damping term explicitly to the desired coefficient of restitution c r (Machado et al., 2012). Starting from the purely elastic HERTZ'ian surface deformation force Kδ n , a non-constant hysteresis damping factor χ is introduced for calculation of the damping part: In order to adapt the damping term to different contact situations and to receive a constant coefficient of restitution, the relative entrance velocityδ (−) of both masses at contact start is additionally used for the calculation of the hysteresis damping factor. Table 3 shows a selection of different definitions of the hysteresis damping factor on the literature (Machado et al., 2012).
In the Modelica implementation the hysteresis damping factor χ can be calculated by a conventional equation, as long asδ (−) is defined at all simulation times including initialisation. However, the initial impact velocitẏ δ (−) has to be assigned within an algorithmic section at discrete time points whenever a new contact starts. Although this model still has discrete and algorithmic parts as well as requires for state events, the numerical integrations were stable, as long as the tolerances were set Oct 08-09, 2020, Tokyo, Japan Table 3. Definition of different hysteresis damping models.

Model Damping factor
FLORES et al.  Table 4). Figure 10 shows the resulting trajectories of contact force versus indentation. The different definitions of the hysteresis damping factors lead to varying energy losses and therefore also varying exit velocities. For the injector, the model of GONTHIER et al. showed the best results in terms of the reached restitution coefficient, compared to the desired one (c r,desired = 0.3, c r,reached. ≈ 0.295).
The specific benefit of all hysteresis damping models is a constant percentage of velocity losses during an impact, which does not vary when used in different impact situations with different masses and initial velocities. Furthermore, during contact the models are neither showing sticking behaviour nor uncontinuous forces. If one mass gets lifted out of a previously resting contact, the consideration of Equation (12) is still necessary, though.

Validation of the dynamic model
In order to validate the accuracy of the complete electromagneto-mechanical Modelica injector model, Figure 11 shows a comparison against a measurement of the Ψ -Icharacteristics of the injector. Due to the lower electric excitation frequency of 7 Hz in the measurement, the eddy current phenomena are not as significant as they are in the real injection cycle. Instead, the missing consideration of static magnetic hysteresis is apparent, however, the influ-  ence of which will get smaller the faster the electric excitation will be. On the other hand, the roundness and symmetric position of the central hysteresis loop indicates an acceptable quality of the combined effects of the total model.

Approaches for needle control
In order to avoid harmful bouncing of armature and valveneedle as well as injecting minimal amounts of fuel in ballistic mode, research on possible alternative control schemes of the injector needs to be done. The principle idea of controlling the trajectory of moving armatures is not new, though. KIRSCHBAUM investigated solenoid intake valves in resonant operating mode and applied dynamic optimisation methods to different types of models, in order to reach the end stop with a predefined velocity (Kirschbaum, 2001). Similar research on conventional intake valves was done by SCHIEDECK, where a heuristic optimisation method provided the parameterisation of a predefined control scheme (Schiedeck, 2003). With a similar aim, two Modelica-specific methods for synthesis of the injector control voltage are given below.
Model inversion Due to the acausality of Modelica, the simplest approach is to invert the actuator model by defin- ing the "output" movement and leaving the necessary voltage supply indefinite, as it is shown in Figure 12, thus making use of the automatic symbolic manipulation of Modelica compilers. The time delaying influence of the electrical (first order) as well as the mechanical system (second order) requires a three times differentiable signal for the armature position. Furthermore, a simplified magnetic model based on look-up tables as well as a modified hysteresis damping model with a constant damping factor is used in order to ease the numerical solution. Figure 13 shows the results of the inverted solution with a ramp for the desired, yet unfiltered armature movement and the corresponding movement s arm for opening the valve. While the acausality of Modelica enables such inverted solutions at low calculation times, the results in case of the examined injector are largely affected by bouncing effects related to the two-stage movement (c.f. Figure 13, t ≈ 0.61 s). Since the inverted model delivers only exact solutions in terms of the specific movement target and does not allow competing objectives, this method is more appropriate for conventional solenoids without body-to-body impacts.
Optimal control Optimisation techniques can handle such problems, as they are formulated as a minimisation of a scalar objective function J, which can account for several conflicting aims. A possible objective function for optimising the trajectory of one full injection cycle can be expressed by a combination of three criteria with their corresponding scaling factors k: In order to guarantee a fully opened cycle as an output of the optimisation, the term J Penalty is used as a penalty, if the armature is not reaching the end stop. The mid term ensures that a certain amount of fuel gets injected; due to the lack of a fluid dynamics model of the injector, this has been approximated by determining a desired time for the valve being opened. The criterium J Impact finally summarises all individual impact velocitiesδ (−) , which take place during the full cycle. In the present work, the complete Modelica injector model has been compiled to a FMU by the use of Open-Modelica and then was imported to JModelica.org. For the optimisation, the control voltage U needs to be discretised over time and used as the variational input. The NELDER-MEAD heuristic optimisation algorithm (Gedda et al., 2012) can be used afterwards for the minimisation of Equation 14, which leads to the results depicted in Figure 14 with the desired smooth armature movement. With further modifications to the objective function, the optimal control problem can be adjusted to nearly every physical demand like, e.g., safely capturing the armature with a predefined first impact velocity. On the other hand the calculation times are considerably higher than the ones of the inverted model, with about 15 h on an 8 core processor 5 . However, many further improvements are imaginable, e.g., the use of a simpler model plus a refinement of the control with the complete model afterwards.

Conclusion
In this paper a new Modelica magnetics library has been presented, which uses an energy-consistent systems modelling analogy. Herewith, the solvability of complex electromagnetic networks including field displacement effects increases in Dymola and OpenModelica with a drawback of higher integration times. Though, future work should concentrate on the specific compilation problems as well as the differences across different compiler versions. In terms of contact mechanics, different models for simulat- ing impacts have been implemented and compared, with the MSL ElastoGap model being less accurate. The hysteresis damping model of GONTHIER et al. showed the best compromise between integration times, ease of parametrisation and physical representation quality. The solenoid model was used for control synthesis, where especially the heuristic optimisation methods are found to be appropriate for the complex engineering demands in mechatronics and especially actuator development. The application and testing of the found control schemes on real solenoids and injectors are left to future work.

Open-source data
The magnetic and contact mechanic libraries are published by the Digital Library Thuringia DBT 6 .