Modeling a Solenoid Driver with Nonlinear Inductive Load Using Circuit Simulation and Magnetic Flux Measurement

Abstract

This paper describes the procedure for creating a electronic simulation model of a solenoid power electronic driver with a nonlinear inductive load. Furthermore, it discusses the electromagnetic interaction between the driver and the load example electromagnetic valve. The consideration of nonlinear effects in the power electronic components MOSFET and diode is particularly important to distinguish their effects from the nonlinear behaviour of the inductive load.

1. Introduction

Different concepts are characterized by electromagnetic systems to apply an electrical current and electrical voltage source to the inductive part of the electromagnetic valve. The interaction between the software-controlled power electronics and the solenoid as a magnetic component is developed based on the requested function with the highest efficiency and accuracy. Refs. [1,2] present typical power electronic circuits with software-controlled function and the results in voltage and current.

There is a series of various automotive concepts to drive electromagnetic injection valves with an electronic power circuit. Ref. [3] describes a pre-energize method, in which the solenoid is electrically energized in advance by taking advantage of attraction. Ref. [4] concerns the development of a high-speed solenoid valve with a 1 ms switching time to control diesel engine electronically. Ref. [5] evaluates simultaneous solutions of the coupled electric, magnetic, and mechanical problems in the dynamic simulation of a solenoid actuator. The represented automotive papers focused on the interaction of the current and the electronic circuit without considering voltage. The evaluation was performed using Finite Element Method (FEM) analysis. This paper is related to both the current and the additional voltage to improve the evaluation result capability. The magnetic flux [6] represents an additional physical information and interpretation of the energy flow and energy conversion between magnetic and power electronic circuit. To analyse and evaluate the impact of the power electronic concept into the magnet system, the simulation can affect different system parameters. To accurately incorporate technical applications into the development process, such as detecting the opening time through electrical control as described in patent [7], the precise modeling of the system and clear differentiation of the effects of individual components are necessary.

In this paper, a circuit simulation model of a test bench is constructed, which can be used to characterize general inductive loads, such as coils or electromagnets. The experimental setup includes the measurement of magnetic flux, and the circuit simulation model was created using LTspice software (v24.0.12) [8]. The analysis is initially focused to geometrically immutable components, disregarding the effects of motion induction and variable magnetic circuits, but can be considered in further steps.

2. Physical Background

According to Faraday’s law, electromagnetic induction is formulated as

$$\begin{array}{c}\hfill {\oint }_{\partial A}\overrightarrow{E}d\overrightarrow{s}=-\frac{d}{dt}{\int }_A\overrightarrow{B}d\overrightarrow{A}=-\frac{d\varphi }{dt}.\end{array}$$

(1)

The line integral over the electric field $\overrightarrow{E}$ yields the induced voltage $u_{ind}$, and the surface integral over the magnetic flux density $\overrightarrow{B}$ yields the magnetic flux $\varphi$ through the conductor loop.

In a coil with multiple windings n, the individual magnetic fluxes $\varphi$ are combined into a linked magnetic flux $\psi =n·\varphi$.

$$\begin{array}{c}\hfill u_{ind}=-\frac{d\psi }{dt}.\end{array}$$

(2)

For a geometrically immutable component (e.g., without armature movement), the magnetic flux depends on both time t and current i, allowing Equation (2) to be rewritten using multivariable chain rule

$$\begin{array}{c}\hfill u_{ind}=-\frac{\partial \psi }{\partial i}\frac{di}{dt}.\end{array}$$

(3)

The negative sign in Equation (3) corresponds to Lenz’s law, leading to the current flowing in the opposite direction to the voltage. If current and voltage are defined in the same direction, the negative sign is eliminated

$$\begin{array}{c}\hfill u_{ind}=\frac{\partial \psi }{\partial i}\frac{di}{dt}.\end{array}$$

(4)

In case of linear behavior, one obtains

$$\begin{array}{c}\hfill u_{ind}=L\frac{di}{dt}.\end{array}$$

(5)

A real inductor consists of the inductive component and the ohmic resistance, as seen in Figure 1.

This leads to

$$\begin{array}{c}\hfill u_{ind}=u-R·i.\end{array}$$

(6)

By separating variables in Equation (2), the magnetic flux can be solved, resulting in the integration of the induced voltage over time

$$\begin{array}{c}\hfill \psi ={\int }_t{u}_{ind}dt={\int }_t(u-R·i)dt.\end{array}$$

(7)

According to Ampère’s law for constant value or low-frequency fields, one obtains

$$\begin{array}{c}\hfill {\oint }_{\partial A}\overrightarrow{H}d\overrightarrow{s}={\int }_A\overrightarrow{J}d\overrightarrow{A}=i,\end{array}$$

(8)

with the magnetic Field $\overrightarrow{H}$, the current density $\overrightarrow{J}$. Equations (1) and (8) allow the formulation of the dependencies

$$\begin{array}{c}\hfill i\sim H,\psi \sim B.\end{array}$$

(9)

Thus, by measuring voltage u and current i, the linked magnetic flux can be determined according to Equation (7). Subsequently, the linked magnetic flux is graphically represented over the current, and conclusions about the magnetic hysteresis can be drawn according to Equation (9).

Equations (5) and (6) combined leads to

$$\begin{array}{c}\hfill u=R·i+L\frac{di}{dt}.\end{array}$$

(10)

The current profile during the switch-on behavior at constant voltage $u_0$ can be calculated as

$$\begin{array}{c}\hfill i\left(t\right)=i_{max}·(1-e^{-\frac{R}{L}t}).\end{array}$$

(11)

After a few rearrangements, the constant inductance L can be determined based on the half-life $t_h$, at which the current reaches $50\%$ of its final value $i_{max}$ ($i\left(t_h\right)=0.5·i_{max}$)

$$\begin{array}{c}\hfill L=\frac{R·t_h}{ln\left(2\right)}.\end{array}$$

(12)

3. Magnetic Flux Measurement

Refs. [9,10] offer an overview of the setup and functionality of fuel injection valve control. The experimental setup used here is based on its application in the automotive sector. In the actual setup, a power supply provides the voltage $u_{batt}$, as can be seen in Figure 2.

A universal solenoid driver (USD) transfers this voltage to the terminals P and N. Between P and N, any inductive component can then be connected. These components can be modeled as a series connection of ohmic resistance and inductance, as mentioned in the previous section. The USD is able to implement arbitrarily controlled current profiles, but in this paper, a simple pulse operation is employed.

The MOSFET $M_1$ controls the current flow i, while the MOSFET $M_2$, on the other hand, is always in the on-state in the application being used here. However, $M_2$ must be considered in the following simulation. One section runs from time $t=[t_0,t_2]$, where $M_1$ is turned on while $\mathsf{\Delta }{t}_1=[t_0,t_1]$, and while $\mathsf{\Delta }{t}_2=[t_1,t_2]$, $M_1$ blocks the current flow from the power supply. By using current clamp A and a differential probe V connected to an oscilloscope, current i and voltage u can be measured from P to N.

The diodes D₁ and D₂ determine the direction of the current flow. From $t_0$ to $t_1$, the current can flow from $u_{batt}$ to ground through $M_1$, $D_1$ and $R_c$ into P, and from N to ground (GND) through M₂, as can be seen in Figure 3. In this case, a magnetic field can build up across the inductive component between P and N, storing magnetic energy within it.

If $M_1$ is switched from the on-state to the off-state at time $t_1$, according to Lenz’s law, a voltage $u_{ind}$ ($u=u_{ind}+u_{ohm}$) is induced between P and N which attempts to maintain the current flow from P to N. Now, the current can flow from GND through $D_2$ and $R_C$ to P, and from N through $M_2$ back to GND.

The last element is the resistance $R_c$. This corresponds to the sum of conductor and contact resistances.

To measure the magnetic flux according to Equation (7), it is necessary to determine the ohmic resistance R of the inductive component. Under the fundamental assumption that the first phase remains activated until the current and voltage of the coil reach saturation at time $t>t^{\prime }$, seen at Figure 4, leads to

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill u& =u_{ohm},\hfill \\ \hfill i& =i_{max},\hfill \end{array}\}\mathrm{for}{u}_{ind}\to 0\mathrm{at}t>t^{\prime }.\end{array}$$

(13)

The time instant $t^{\prime }$ is chosen such that the changes in voltage and current approach zero. The closer $t^{\prime }$ is set towards $t_0$, the less accurately R is determined, leading to errors in the $\psi$-i characteristics. This allows the ohmic resistance R to be calculated using Ohm’s law

$$\begin{array}{c}\hfill R=\frac{u}{i}(t>t^{\prime })=\frac{u_{ohm}}{i_{max}}.\end{array}$$

(14)

The advantage of this method is that the resistance can be automatically determined for each cycle k ($t_0\to t_2$) during the measurement. This results in very good repeatability, as measurement inaccuracies and temperature drift of R are immediately taken into account. The disadvantage is that current and voltage must go into saturation.

A further method is described in patent [11]. In this approach, a parameter variation is carried out for each cycle k, and for various values of $R_k$, the $\psi$-i characteristics are determined. The value of $R_k$ that results in a closed $\psi$-i characteristics is then the ohmic resistance. Over multiple cycles k, the value can then be determined as

$$\begin{array}{c}\hfill R=\frac{1}{k}\sum R_k.\end{array}$$

(15)

This allows for the determination of resistance even when current and voltage do not reach saturation. However, in this case, R is determined before the actual measurement of $\psi$ through a series of measurements.

Due to its repeatability, the first method for resistance determination is employed in subsequent measurements.

A solenoid valve with a fixed armature is used as the test inductance between the terminals P and N. The components of the solenoid valve that participate in the electromagnetic force effect are shown in Figure 5.

The coil carrier does not contain any magnetizable materials, so it can be viewed as an air coil without a built-in core. On the other hand, the anchor and core, which are made of the same material, can be magnetized and show hysteresis behavior. To characterize the material properties of the armature and core, the armature is fixed in such a way that it closes the magnetic circuit with the core and the air gap becomes zero. The fixation is completed using spacer plates that are installed on the valve side, as can be seen in Figure 6.

To model the diodes and MOSFETs from Figure 2, the coil is measured without core and armature. From this, a constant value for the inductance and ohmic resistance of the air coil can be determined. This avoids the influence of nonlinear inductance, making modeling of the semiconductor components easier. The value of the inductance for the air coil L can be determined by measuring the half-life and calculated according to Equation (12). After the semiconductor components have been modeled, the armature and core can be inserted into the coil and their nonlinear properties can be examined. It is important to note that the sudden switching on of the voltage creates eddy currents in the armature and core, which in turn generate a counter-field. Therefore, the resulting $\mathsf{\Psi }$-i characteristic curve includes not only nonlinear hysteresis properties but also eddy current effects. However, since only an $RL$ series circuit is being used, and the inductance includes all nonlinear effects, this circumstance can be neglected.

The complete experimental setup is shown in Figure 7.

4. Magnetic Flux Simulation—USD-Modelling

Next, a simulation model is developed for the measurement, a kind of a digital twin of the real experimental setup. In context of this paper, a circuit simulation of the measurement setup from Figure 2 was implemented in LTspice. First, the results of the final equivalent circuit are presented, followed by a description of the modeling of the semiconductor components. The voltage $u_{batt}$ is increased from 3 V to 24 V in steps of 3 V. The comparison between the measured and simulated results in current i and voltage u, as illustrated in Figure 8, shows a good alignment.

First, the modeling of the diodes is considered. These can, for example, be modeled using two different approaches. A simple approach involves a piecewise linear approximation of the diode characteristic curve with the diode voltage $u_D$. Below the threshold voltage $u_{th}$ ($0<u_D<u_t$), no voltage drop occurs across the diode ($u_D=0$). Above the threshold voltage ($u_D>u_{th}$), an ohmic value is defined that linearly increases the voltage drop across the diode. According to the datasheet, the threshold voltage is $u_{th}=0.4$ V. A more accurate approach would be a nonlinear modeling as shown in [12,13]. In the range of diode voltage $u_D>0$, where they are primarily operated, the nonlinear behavior is described by the conductivity depending on the diode voltage, as can be seen in Equation (16).

$$\begin{array}{cc}\hfill G& =\left(\frac{i_s}{n·u_t}\right)·exp\left(\frac{u_D}{n·u_t}\right),\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill u_t& =\frac{k_B·T}{e},\hfill \end{array}$$

(17)

with $i_s$ as saturation current, n as the emission coefficient and $u_t$ as thermal voltage depending on Temperature T, calculated with Boltzmann constant $k_B$ and elementary charge e. Lastly, the zero-bias junction capacitance $Cjo$ was included in the diode-modeling. With the exception of the parameter n, all other values can be extracted directly from the datasheet. The emission coefficient is then chosen such that the resulting diode characteristic curve matches the curves from the datasheet. The datasheet provides a maximum value for the saturation current. Therefore, in addition to the maximum value, a reduced value was also tested. The latter provides a better approximation of the characteristics, as well as a better agreement between the simulation results and the measurement results. Hence, the reduced value of $i_s$ was finally adopted into the model. In Figure 9, one can observe the individual approaches in a direct comparison of the diode characteristic curves.

In Figure 10, the current and voltage profiles are presented, showing the impact of various diode models. Notably, the voltage profile during freewheeling exhibits a significant influence.

The MOSFETs were modeled based on a VDMOS model extracted from the LTspice library, resembling a similar MOSFET. Subsequently, individual parameters were adjusted according to the datasheet based on various calculations found in [13]. In the saturation region of the MOSFET, the drain current $i_D$ is governed by the equation

$$\begin{array}{cc}\hfill i_D^{sat}& =\frac{\beta }{2}(u_{GS}-u_{th}{)}^2,\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill \beta & ={\mu }_e{C}_{ox}\frac{W}{L}.\hfill \end{array}$$

(19)

In Equation (18), $u_{GS}$ represents the voltage across the gate and source, and $u_{th}$ is the threshold voltage. $\beta$ is a component-specific parameter and, according to Equation (19), depends on the charge carrier effective mobility ${\mu }_e$, gate length L and width W, and the gate oxide capacitance per unit area $C_{ox}$. Based on the $i_D-u_{DS}$ characteristic from the datasheet with various curves for $u_{GS}$, two points in the saturation region can be selected and substituted in to Equation (18)

$$\begin{array}{cc}\hfill i_{D1}^{sat}& =\frac{\beta }{2}(u_{GS1}-u_{th}{)}^2,\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill i_{D2}^{sat}& =\frac{\beta }{2}(u_{GS2}-u_{th}{)}^2.\hfill \end{array}$$

(21)

Dividing Equation (20) by Equation (21) and solving for $u_{th}$ allows determination of the threshold voltage

$$\begin{array}{c}\hfill u_{th}=\frac{u_{GS1}\sqrt{i_{D2}}-u_{GS2}\sqrt{i_{D1}}}{\sqrt{i_{D2}}-\sqrt{i_{D1}}}.\end{array}$$

(22)

The static Drain-Source on-resistance $R_{on}$ can be directly extracted from the datasheet, while all other resistances were adopted from the similar-resistance in the LTspice library. The switching characteristics of MOSFETs exhibit a strong dependence on the capacitances of the MOSFET. The capacitances required for the LTspice-model can be determined using the input capacitance $C_{iss}$, output capacitance $C_{oss}$ and reverse transfer capacitance $C_{rss}$ specified in the datasheet. The capacitance between Gate and Source $C_{GS}$, Drain and Source $C_{DS}$ and between Drain and Gate $C_{DG}$ can be determined by

$$\begin{array}{cc}\hfill C_{GS}& =C_{iss}-C_{rss},\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill C_{DS}& =C_{oss}-C_{rss},\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill C_{DG}& =C_{rss}.\hfill \end{array}$$

(25)

Finally, the transconductance parameter $g_m$ was determined based on the slope of the $i_D-u_{GS}$ characteristic curve by

$$\begin{array}{c}\hfill g_m=\frac{\partial i_D}{\partial u_{GS}}.\end{array}$$

(26)

However, there remains a residual deviation, as the unknown parameters, not included in the datasheet, sourced from a similar MOSFET in the LTspice library.

After determining the semiconductor components, the only remaining parameter is the contact resistance. This resistance is significantly higher in value than the ohmic resistances of the diodes and MOSFETs in the forward direction. This can be relatively easily determined using a voltage divider, by connecting a resistance at terminals P and N, both in simulation and measurement. The value of $R_c$ only needs to be adjusted so that the voltage drop across P and N matches the measurements.

5. Magnetic Flux Simulation—Nonlinear Inductivity

If the components of the USD circuit are described with sufficient accuracy, the air coil as a linear inductance can be replaced with a nonlinear component, in this case, the valve from Figure 5 and Figure 6. By varying the voltage $u_{batt}$ of the power supply from 3 V to 24 V in 3 V increments, the $\mathsf{\Psi }$-i characteristic curve is determined for each voltage-step. As the core and armature are constructed from ferromagnetic material, they exhibit hysteresis loops in the B-H characteristic curve. Therefore, according to the dependency (9) for the $\mathsf{\Psi }$-i characteristic curve, a nonlinear hysteretic course arises. The enclosed area represents demagnetization losses. As mentioned earlier, due to the pulsed switching behavior, eddy currents occur in the conductive core and armature, further broadening the enclosed area. The varying voltages consequently result in different currents, which, according to Equation (8), are the cause of the magnetic field. Both the strength and the pattern of the magnetic field are influenced by the material properties. As seen in Figure 11, the flux curves exhibit minimal curvature at lower voltages. In these regions, the hysteresis loop can often be approximated with constant inductances. Raising the voltage pulses further reveals the typical nonlinear curvature, extending into the saturation region.

Now, a nonlinear $\mathsf{\Psi }$-i characteristic curve is to be found to approximate the measured characteristics for all defined voltages. For this purpose, the respective flux maxima are used as nodes for interpolation including the starting point $(0,0)$. Polynomial functions are used as the basis for interpolation, once of the third order and once of eighth order. The coefficient for the constant term is set to 0 to ensure that the starting point is perfectly matched. The interpolation was performed using the MATLAB^® curve-fitting toolbox (version R2024a). Figure 12 displays the resulting curves, along with the deviation from the nodes.

The eighth-order polynomial approach generally provides a better fit to the nodes than the third-order. It is important to notice that the polynomial functions are only valid in the range from 0 V to 24 V. Beyond this range, they exhibit a significantly different curve, leading to unphysical results. The comparison between the measured and simulated results with 8th-order polynomial approach shows a good alignment, as can be seen in Figure 13.

In the following figures, the course of the magnetic flux for two voltage-nodes is detailed. Figure 14 demonstrates a very good agreement of the eighth-order with the measurement, while Figure 15 serves as an example where the deviation is greater.

The method using flux maxima as nodes for interpolation performs excellently for the rising flank of the hysteresis and adequately in the steady state. However, deviations are noticeable in the falling flank. This is indicated by the fact that while the eighth-order provides a better approximation according to Figure 12, it consistently exhibits a poorer performance in the falling flank compared to the approach using the third-order. This can be attributed to the slightly different curvature behaviors of the two approaches. While the eighth-order better corresponds to the curvature of the rising hysteresis flank and captures the steady-state endpoint more accurately, the third-order exhibits better alignment with the curvature behavior of the falling hysteresis flank.

6. Conclusions

In summary, an electric system for controlling inductances was modeled as a circuit simulation. To accomplish this, the system was initially loaded with a constant inductance to avoid nonlinear effects from conductive, ferromagnetic materials, while neglecting motion induction. As a result, the nonlinear effects of the semiconductor components could be more easily examined and modeled with greater precision. During this process, the simulation results were compared with measurements, and parameters were adjusted accordingly. Once the semiconductor components were accurately modeled, the constant inductance was replaced by a nonlinear, hysteresis-dependent inductance. By gradually increasing the voltage, multiple flux characteristics were recorded, and through interpolation using polynomial approximations of various orders, nonlinear functions for the magnetic flux was determined. These polynomial functions, which are nonlinear but independent of hysteresis, are intended to approximate the measured flux characteristics dependent on hysteresis.

Since a hysteresis-independent approach was used for modeling the inductances, this approach can only be applied to soft magnetic materials with low demagnetization losses. For hard magnetic materials, only either the rising or falling hysteresis flank can be approximated accurately. Furthermore, eddy current effects were neglected. To account for these in the circuit simulation, a parallel resistance to the inductance can be used to represent eddy current losses [14]. Also, as noted at the beginning, motion induction was set to zero due to the fixed anchor. The influence of motion induction leads to

$$\begin{array}{c}\hfill \frac{d\mathsf{\Psi }\left(i\right(t),x(t\left)\right)}{dt}=\frac{\partial \mathsf{\Psi }}{\partial i}\frac{di}{dt}+\frac{\partial \mathsf{\Psi }}{\partial x}\frac{dx}{dt}=\frac{\partial \mathsf{\Psi }}{\partial i}\frac{di}{dt}+\frac{\partial \mathsf{\Psi }}{\partial x}{v}_x\left(t\right).\end{array}$$

(27)

The term appearing in Equation (27) for the motion induction requires further modeling. Firstly, an approach for the magnetic flux needs to be found, which accurately describes the dependence on displacement in addition to the current. Additionally, the resulting velocity needs to be determined through a coupled mechanical model, depending on the occurring magnetic force. If a more detailed description of hysteresis including minor loops is necessary, a more comprehensive model like the Preisach model must be employed, described in [15,16]. In [17], various hysteresis models are compared with each other in terms of error, complexity of parameter estimation, and the effort required for validation through measurement. The respective strengths and weaknesses of the models are presented.

A simpler solution approach to consider eddy current effects and motion induction is offered by transferring the problem into a Finite Element Method (FEM) model. For this purpose, in the circuit model shown in Figure 2, no resistive and inductive elements are defined as loads. Instead, the terminals P and N are coupled with a 2D or 3D FEM model. The FEM model itself must accurately represent the real geometry, and instead of $\mathsf{\Psi }$-i characteristics, B-H characteristics are defined as material properties for the core and armature.