An Alternative Multi-Physics-Based Methodology for Strongly Coupled Electro-Magneto-Mechanical Problems

Abstract

The analysis of complex systems tends to be approached through a separation and a simplification of the main macro phenomena and, thus, addressed through dedicated techniques, tools, and algorithms. A smart and interesting possibility, instead, is represented by the so-called model-based design analysis, which allows one to interface phenomena coming from interactions of different physical natures. This paper aims to propose a multi-physics Matlab/Simulink^®-based architecture that allows one to integrate general and strongly non-linear coupling phenomena, taking efforts from two novel implemented bi-directional co-simulation routines based on Spice^® and ESRF Radia^® engines. Emphasis is dedicated to the discussion and description of the co-simulation algorithms and processes characteristic of these routines, which allow the analog electronic and the magneto dynamic domain’s integration under a single simulation environment. To highlight the reliability of the multi-domain architecture and to validate the reported co-simulation results, a comparison with the experimental measures obtained on an innovative MEMS electromagnetic actuator are proposed.

1. Introduction

Time-domain analysis represents a prerogative for the comprehension and design of mechatronic devices. Often, the complexity of the numerical modeling of such systems is narrowed and simplified through the isolation of the main macro-phenomena, and thus is solved through dedicated languages, tools, and algorithms [1,2]. However, in applications where the strict interdependence between the different physical domains and scales cannot be decoupled, such as mechatronic devices, strongly coupled electro-magnetic applications (WPT, electric vehicle re-charging, etc.), and active analog electronics controllers, the introduction of a complete multidomain co-simulation architecture is required [3,4,5,6,7].

From a holistic and practical point of view, one of the main difficulties arises from the integration of the different platforms on which the partial models are developed. Among the others, Matlab/Simulink^® represents one of the most advanced graphical programming environments for the simulation of dynamic systems, being developed for the analysis of multipurpose and general nature problems [8]. In this study, we propose a novel Matlab/Simulink-based architecture for multi-physics analysis, devoting particular emphasis to the interaction between the mechanical, magnetic, and analog electronic domains, taking efforts from two novel implemented bi-directional co-simulation routines based on Spice^® and ESRF Radia^® engines. In detail, this paper describes the logics and the sub-processes executed by the novel co-simulation algorithms for bi-directional communication of the effort and flow variables incoming and outgoing from the Matlab/Simulink leading routine and the analog electronic, and the magneto-static/dynamic ambient solved through equivalent circuit analysis (ECA) and 3D finite volume modeling (3DFVM), respectively.

To validate the proposed methodology, an innovative MEMS-based electromagnetic actuator was taken under analysis, and the numerical results obtained by the multi-physics model were compared with the experimental relatives reported by Lee et al. in [4]. In the following, Section 2 introduces the whole co-simulation architecture, highlighting the main interexchange processes and the related timing characteristics of the different domains. Importance is given to the two novel bi-directional co-simulation algorithms for the integration of the analog electronic and magneto-static domain through PySpice-Simulink and ESRF Radia-Simulnk routines, faced in Section 3 and Section 4, respectively. Section 5 highlights the potentiality of the proposed multi-domain architecture through a validation process based on the experimental results obtained from an electromagnetic actuated MEMS proposed in the literature. Finally, appropriate conclusions are drawn in Section 6.

2. The Proposed Multi-Domain Architecture

Figure 1 depicts a block diagram representation of the proposed implementation logic behind the co-simulation multi-scale and multi-physics modeling technique. The core of the model, developed in Matlab/Simulink, discretizes the time of the simulation interfacing the different physical domains.

This discretization process represents one of the key aspects of the co-simulation algorithm. It takes place independently on four different timelines. In detail, as shown in the lower part of Figure 1, each physical domain has its own simulation time window and a characteristic sampling frequency. The Matlab/Simulink leading ambient defines the maximum temporal length for the entire simulation, discretizing it through a variable step integrator. The length of this step represents the maximum length of the instantaneous simulation window of the individual physical domains, which are called independently according to their respective calling frequencies; these are necessarily lower than or equal to the minimum leading one. In turn, as in the case of the electronic domain call, the Spice integrator further discretizes the relative timeline, with frequencies in the order of ${10}^9$ Hz, returning at the end of the co-simulation process, for the considered Matlab/Simulink time instant the results obtained from the last Spice discretization interval. This aspect, firstly, allows independent integration of the three standard tools for dedicated analysis as well as to sample/call the different domains at an optimal frequency, thus optimizing the computational effort of the CPU and preventing the aliasing phenomenon.

Inside the MATLAB/Simulink master routine, two other main functions are also implemented: the high-frequency parameter estimator, which computes the induced variation in parameters and variables affected by the interaction with different physical domains, and the firmware and control logic routine, which allows one to embed a firmware control loop in the case that the studied device is governed by a microcontroller. In detail, a complete multibody model of the full device assembly under analysis is implemented via MSC Adams^® and co-simulated through Adams control, the Mathematica open-source package: ESRF Radia solves the static and dynamic electro-magnetic fields through a 3D FVM algorithm, and finally the equivalent circuit properties are integrated and solved through the Spice engine. The whole dynamics behavior is integrated at every time step and updated; the resulting sets of PDE/ODE of each domain are settled, taking effort off the dedicated platform solvers; and finally the bi-directional communication is taken in charge by two novel implemented and dedicated Matlab/Python routines (see Section 3 and Section 4).

3. Spice-Based Analog Electronic Domain Integration

3.1. Introduction

The increasing demand for accuracy and reliability in the design and analysis of complex electronic integrated systems has, in recent years, pushed toward a greater use of device-level circuit simulators. Their strengths are based on the conformity to fundamental physical principles [6], which in turn introduces long-time simulation, losses in computing performances, and difficulty in multi-domain integration. On the contrary, a system-level ambient such as Matlab/Simulink can guarantee high computational performances and multipurpose analysis at the expense of a reduced range of possible simulations, due to lack of transient, AC/DC, or sweeps analysis capabilities and a limited availability of complete components libraries [9].

The proposed work therefore introduces and defines a new co-simulation methodology between a device-level circuit simulator (Spice) and a system-level ambient (Matlab/Simulink) so as to ensure high fidelity in the simulation of the related physical phenomena and to the broad range of events that can be studied.

In the literature, few attempts have been made toward this direction; among the others, a work was proposed by Madbouly et al. in [10] in which a set of simulation tools called MATSPICE were presented that introduced a preliminary step for the optimization of parameters and the rapid iteration of circuit resolutions through a SPICE-based solver engine. Another significant implementation was proposed in [8], where authors described an innovative and interesting approach to the problem through a Simulink-based S-function implemented as a Level-2 M-file capable of coupling an NGSpice engine with a Matlab/Simulink model.

The proposed work instead aims to integrate a device-level simulator such as a Spice-based circuit at the system-level, such as a Matlab/Simulink environment, through a completely different approach. A dedicated Simulink/Matlab-based function was executed through the Matlab Engine API for Python, the implemented Python-based routine, which integrated the open source module PySpice [11], which in turn allowed the NGSpice/Xyce simulator to be run in server mode, as well as the bi-directional communication of the effort and flows variables among the co-simulated ambientes.

3.2. The Proposed Spice-Based Co-Simulation Algorithm

Figure 2 depicts the logic behind the algorithm of co-simulation between the Spice-based circuit integrator and Matlab/Simulink.

The flow chart reports two main macro steps: the initialization process and the iteration cycle. The circuit initialization starts with a fast implementation of the circuit schematics through a Graphic User Interface (GUI) editor. The components network is grouped inside the so called “.net” file NETLIST. To be correctly integrated and executed by the co-simulation algorithm, the NETLIST needs to be modified according to what is reported in Figure 3.

The title section marked in the picture with 1 was set to the first line since this is never executed by the Spice solver. Marker 2 highlights the core of the NETLIST, where all the components network is defined, based on the circuit nodes’ connections. As reported, the constitutive values of some components are defined in brackets; these bracketed parameters represent the turning point of the co-simulation process. In fact, as listed by the section numbered with marker 3, the so called “.param par…” defines the INPUT values that are passed from Simulink to Spice during the co-simulation exchange and allow a continuous modification of the circuit caused by the interaction with the other physical domains. This step would be extremely complex or even not replicable through a purely device-level-based approach.

As previously described and shown in Figure 1, during each Simulink call, the Spice integrator generates a new independent timeline that needs to be temporally relocated in the main macro timeline and initialized on the basis of the previously obtained results. These initial condition parameterizations are listed in the section marked as 4 through the so-called “.param ic…” and “.ic…”. Finally, Section 5 defines and lists all the equivalent NGSpice/Xyce-based models of the real electronic components.

The so-initialized circuit model is ready to be co-simulated on the basis of the following steps: (I) Simulink calls off the communication routine framework through the Matlab Engine API for Python that executes the Spice solver based on the open source PySpice module; (II) NGSpice/Xyce solves the circuit NETLIST; (III) the routine extracts the last V/I Spice discretization results, as Spice samples the electronic phenomena at a higher frequency (i.e., 1 GHz) compared with the Simulink leading routine call; (IV) a Matlab function re-initializes the circuit NETLIST updated on the last extracted Spice results and on the solved ones and updates INPUT parameters coming from the interaction with the other domains (e.g., mechanical, magnetic…); and (V) the loop is closed, transferring the resulting effort and flow variables to the co-simulated domains.

4. ESRF Radia-Based Magnetostatic Domain Integration

4.1. Introduction

The other key structure developed for the presented electro-magneto-mechanical architecture is constituted by the ESRF Radia-based bi-directional communication algorithm. In the literature, there are several interesting case studies reported, mainly focused on the field of motor drives or magneto-strictive actuators. It is important to remark how these reported works face the analysis of the electro-magnetic interaction through a finite element analysis (FEA) approach [12,13].

In our study, we propose an alternative perspective with respect to this classic FEM-based analysis, taking efforts to consider the potentiality of the finite volume modeling (FVM) [14].

4.2. The Proposed ESRF Radia-Based Co-Simulation Algorithm

As for the co-simulation algorithm reported in Section 3 also for the case of magneto-static integration, the whole procedure is characterized by the same two main macro steps: the initialization process and the iteration cycle.

The 3D Finite Volume Model (3DFVM) was first implemented in Wolfram Mathematica^® language and then imported into a Python module. The script now starts its initialization process to be suitable for the co-simulation analysis. Now that the model is fully set, it is passed to the iterative cycle according to what is reported in Figure 4. The main algorithm framework uses a logic that exploits a closed-loop process based on the following steps: (I) the leading Simulink routine calls off and updates, through the Matlab Engine API for Python, the 3D FVM model implemented in Mathematica language. This is integrates into a Py-module and interpreted by the Wolfram Client Library for Python; (II) the ESRF Radia engine solves the model based onto the achieved configuration; (III) the routine extracts the obtained results (e.g., magnetic fields intensity, magnetic forces/torques, etc.); and finally (IV) the loop is closed, transferring the resulting magnetic effort and flow variables to the other physical domains.

From the simulation point of view, the magnetostatic analysis can be processed through two different main strategies: the real-time co-simulation or the pre-simulation, and look-up table interpolation.

In the case that the real application of the device presents a strict correlation between the coupling aspects with two or more kinematic variables, a real-time co-simulation is considered necessary. On the contrary, if the dependency presents a single kinematic correlation, it is possible to pre-simulate the whole 3D magnetostatic behavior for all possible 1Ds configurations and to pass the resulting trends to the co-simulation model in the form of an interpolated curve. This last scenario allows one to considerably increase the computation performances of the entire multi-domain model.

5. The Case Study

5.1. Introduction

To fully complete the description of the main capabilities of the proposed co-simulation architecture and to highlight the potentiality as well as the reliability in a real case scenario, a simple application case study together with a detailed experimental comparison is presented. The experimental results taken into analysis for the numerical validation were extracted from the work reported by Lee et al. in [15], based on an innovative MEMS electromagnetic actuator.

The authors in their work proposed a study about the design, the fabrication, and characterization of a novel micro electro-magnetic actuator comprising a PDMS diaphragm, a polyimide-coated copper micro coil, and a permanent magnet.

5.2. The Device

Ref. [15] analyzed three different devices based on different geometries and characteristics; the proposed experimental-numerical comparisons were based on the results obtained from the device that developed the most effective performances, the so-called Type I design.

In detail, the studied device (see Figure 5) was structured as follows: a $20\mathsf{\mu }\mathrm{m}$ thick copper planar micro-coil was mounted over a glass-polyimide substrate. The coil presented an internal diameter of $2500\mathsf{\mu }\mathrm{m}$ and an external one of $5087\mathsf{\mu }\mathrm{m}$, the in-planar spacing and track width were $80\mathsf{\mu }\mathrm{m}$ and $100\mathsf{\mu }\mathrm{m},\mathrm{respectively}$. A PMMA plate covered the coil-substrate sandwich; at the center of the plate, a circular hole was obtained with a diameter of $4\mathrm{mm}$ that constituted the device’s pumping chamber. On the plate’s upper surface was glued the PDMS diaphragm–magnet coupled assembly. This assembly was constituted by a 1.5 × $1.5\mathrm{mm}$ cylindrical NdFeB permanent magnet glued over a $100\mathsf{\mu }\mathrm{m}$ thick circular PDMS diaphragm with a diameter of $2000\mathsf{\mu }\mathrm{m}$. Further and deeper constitutive information can be found in the reported article [15].

5.3. The Governing Equations

The implementation of a complete numerical model requires one to describe the overall coupled phenomena. The main physical domains for the treatment of the device under analysis are the mechanical domain, the magnetic domain, and the electronic domain. The co-simulation architecture described in the previous sections is here integrated with a set of governing equations to provide a complete and exhaustive modeling framework. This section encompasses the main purposes of this article, i.e., the integration of complex co-simulation tools in order to implement and solve a complete analytical/numerical multi-domain model.

5.3.1. The Planar Micro-Coil Characterization

The electrical DC resistance of the copper micro-coil was defined by authors as equal to $2\mathsf{\Omega }$, due to the field of application and the necessity to take into account skin and proximity effects; therefore, the actual parameterization value of the electrical resistance was identified as the equivalent resistance value for AC excitations.

Shen et al. in [16] proposed a detailed mathematical modeling of the Dowell’s derived equations of an equivalent resistor to introduce in circuit simulation.

The equivalent AC resistance expression $R_{AC}$ was formulated as a combination of different hyperbolic functions $\zeta$ and $\xi$.

$$R_{AC}=R_{DC}{F}_r=R_{DC}\mathsf{\Delta }\left[\mathsf{\zeta }+\frac{2}{3}\left(p^2-1\right)\xi \right]$$

(1)

The proximity $\xi$ and skin $\zeta$ effect [17] instead need to be numerically calculated according to Equation (2), through a combination of harmonic functions of the penetration ratio $\mathsf{\Delta }$.

$$\xi =\frac{\sinh \mathsf{\Delta }-\sin \mathsf{\Delta }}{\cosh \mathsf{\Delta }+\cos \mathsf{\Delta }}\zeta =\frac{\sinh \left(2\mathsf{\Delta }\right)+\sin \left(2\mathsf{\Delta }\right)}{\cosh \left(2\mathsf{\Delta }\right)-\cos \left(2\mathsf{\Delta }\right)}$$

(2)

According to the set of Equations reported in Equation (3), the penetration depth $\mathsf{\Delta }$ is essentially a function of the physical and geometric parameters of the material.

$$\mathsf{\Delta }=\frac{\sqrt{\eta }{d}_w}{\delta }\eta =\frac{t\sqrt{k}{d}_w}{h_c}\delta =\sqrt{\frac{\rho }{\pi \mu f}}$$

(3)

where

• $\eta$ is the porosity factor of the material.
• $d_w$ is the equivalent thickness of the material and for a round wire and is equal to $\frac{\sqrt{\pi }{d}_i}{2}$, where $d_i$ is the diameter of the wire.
• $\delta$ is the skin dept in the wire.
• $t$ is the number of turns per layer.
• $k$ is the number of Litz wire, which is 1 for round wire.
• $h_c$ is the height of the core window.
• $\rho$ is the resistivity of the wire material.
• $\mu$ is the permeability constant equal to ${\mu }_r{\mu }_0$.
• $f$ is the principal operating frequency.
• $p_i$ is the equivalent number of layers.

A more in-depth description and identification of the parameters as well as the whole mathematical demonstration can be found in [16].

The other two main parameters needed for the characterization of the micro-coil were the stray capacitance and inductance. Costa proposed in his work [18] a complete description of the parasitic capacitance behavior of a planar micro-coil. Starting from the geometry’s properties, such as the number of turns $n$, the internal $R_{in}$, external radius $R_{out}$, and the distance between two turns $d=\left(R_{out}-R_{in}\right)/n$, the author proposed a set of equations for the analytical formulation of the equivalent stray capacitance. Modeling the whole planar coil as a series of i-th parasitic capacitances delimited by the polar coordinates $a$ and $b$ defined in Equation (4) with the potential and the electrical flux between two turns, it is thus possible to define the global $X$ and $Y$ constitutive Equations (5) and (6).

$$a=R_{in}+\left(i-1\right)d;b=a+id;fori=\left\{1\le i\le n-1\right\}$$

(4)

$$X=\frac{-b+\ln \left(b^{3}a\right)}{2a^2}+\frac{1}{4a^2}+\frac{1}{2b}-\frac{\ln \left(b^2\right)}{b^2}-\frac{1}{4b^2}$$

(5)

$$Y=-3\ln \left(\frac{b}{a}\right)+b$$

(6)

The effective height of the coil $z$ is defined in Equation (7).

$$z=8.25\times {10}^3\left(n-10\right)+0.1401$$

(7)

Plugging the set of reported equations, it is possible to formulate the equivalent stray capacitance $C$ as the sum of the i-th coupled capacitors $C_i^{\prime }$.

$$C^{\prime }=\frac{4z\pi ϵ}{a^2}\frac{Y}{X}$$

(8)

$$C={\displaystyle \sum }_{i=1}^{n-1}{C}_i^{\prime }$$

(9)

The inductive behavior instead is modeled through the mathematical formulation proposed by Mohan et al. in [19]. Herein are proposed different methods for the analysis of the planar coil. The method adopted during the research activity was the so-called Current Sheet Expression. Where $c_1,c_2,c_3$ and $c_4$ are experimentally tabled parameters, $d_{avg}=\left(d_{out}-d_{in}\right)/2$ is the average diameter of the coil, $n$ represents the number of turns, ${\mu }_0$ is the magnetic permeability of free space, and $\rho$ is the so-called fill ratio.

$$\rho =\frac{d_{out}-d_{in}}{d_{out}+d_{in}}$$

(10)

$$L=\frac{\mu n^2{d}_{avg}{c}_1}{2}\left(\ln \left(\frac{c_2}{\rho }\right)+c_3\rho +c_4{\rho }^2\right)$$

(11)

The resulting R-L-C parameter estimations are listed in

Table 1. Estimated R-L-C parameters @ 100 Hz.

Parameter	Value	Unit
AC resistance: $R_{AC}$	2	$\left[\mathsf{\Omega }\right]$
Stray capacitance: $C$	$1.36\times {10}^{-10}$	$\left[\mathrm{F}\right]$
Stray inductance: $L$	$4.766\times {10}^{-7}$	$\left[\mathrm{H}\right]$

.

5.3.2. The PDMS Diaphragm Characterization

The other important characterization process concerns the analysis of the dynamic mechanical behavior of the PDMS membrane properties. Figure 6 depicts an electro-magnetic actuated mass-spring-damper equivalent schematic that represents the model adopted for the characterization of the main dynamic parameters. The here-reported model can be fully defined based on three main lumped parameters: the equivalent spring $k$, the equivalent damping $c$, and the equivalent mass $m$.

The total equivalent mass was obtained from the CAD model of the device and was the result of the sum of the PDMS diaphragm mass plus the mass of the NdFeB permanent magnet. The lumped value of the equivalent spring was instead obtained from the reduced static equilibrium at the mass of the system Equation (12).

$$F_{mag}=F_g+F_{el}$$

(12)

where $F_g$ is the weight force of the total mass, $F_{el}$ is the force exerted by the equivalent spring, and $F_{mag}$ is the Lorentz formulation of the magnetic force exerted by an n-turning coil. The resulting expression of Equation (12) is reported below in Equation (13).

$$nli\mathsf{\Lambda }B=\left(m_{diaphragm}+m_{magnet}\right)g+k{w}_{max}$$

(13)

Coupling the resulting formulation with the theory of diaphragm deflection (see Equation (14)) proposed by Bao in [20], it is possible to trace the formulation of the equivalent winding length and the equivalent stiffness parameter, reported in the set of Equation (15).

$$w_{max}=\frac{F{r}^4}{S64D}D=\frac{E{h}^3}{12\left(1-{\nu }^2\right)}$$

(14)

$$l=\frac{w_{max}\pi r^{2}64D}{NiB{r}^4}k=\frac{nliB-\left(m_{diaphragm}+m_{magnet}\right)g}{w_{max}}$$

(15)

The constitutive parameters’ nomenclature and the relative values are listed in

Table 2. Measured and estimated parameters.

Parameter	Value	Unit
Diaphragm radius: $r$	1000	$\left[\mathsf{\mu }\mathrm{m}\right]$
Diaphragm thickness: $h$	100	$\left[\mathsf{\mu }\mathrm{m}\right]$
Diaphragm surface: $S$	3.15	$\left[{\mathrm{mm}}^2\right]$
PDMS Young’s modulus: $E$	550	$\left[\mathrm{kPa}\right]$
PDMS Poisson’s ratio: $\nu$	0.5
Coil current: $i$	0.6	$\left[\mathrm{A}\right]$
Number of turns: $n$	10
Magnetic flux density: $B$	0.0416	$\left[\mathrm{T}\right]$
Max deflection @ $0.6A$: $w_{max}$	150	$\left[\mathsf{\mu }\mathrm{m}\right]$
Equivalent coil length: $l$	7384	$\left[\mathsf{\mu }\mathrm{m}\right]$
Equivalent spring stiffness: $k$	0.0078	$\left[\mathrm{N}/\mathrm{mm}\right]$

.

Due to the lack of information provided in the referenced article [4], the equivalent damping $c$ of the device was tuned numerically and obtained a value equal to $5.53\times {10}^{-6\frac{\left(\mathrm{Ns}\right)}{\mathrm{mm}}}.$

5.4. The Multi-Domain Model

Figure 7 depicts the co-simulation architecture, in which the main effort and flow variables exchanged according to the logic described in Section 2 are highlighted and reported.

The leading routine discretizes the time of the simulation, managing the inter-exchange of variables among the different domains as well as estimating the parameters needed in the analysis based on the set of equations described in the previous paragraphs. In detail, at every time step, the AC resistance $R_{AC}$, the equivalent inductance $L$, and the stray capacitance $C$, together with the voltage amplitude $V$ of the generator, the voltage amplitude of the equivalent back e.m.f. generator $V_{backe.m.f.}$, and the amplitude of the equivalent microcontroller controlling the MOSFET $V_{gate}$ are passed to the analog electronic domain that solves the instantaneously configured circuit and outputs the voltages’ drop across a shunt resistor $\left(V_{n5}-V_{n6}\right)$. According to Ohm’s law, this voltage drop is used to compute the current circulating in the coil. The MBD model outputs the deflection amplitude $w$ and the related membrane deflection velocity $dw$. The deflection amplitude measure is passed to the magnetic domain to update the equivalent magneto-static model, which solves the 3D FVM analysis and passes back the estimated magnetic field density $B$ to the leading routine for the computation of the actuating magnetic force $F_{mag}$. This computed force value is passed back to the MBD system for the dynamic actuation of the device. The deflection velocity, instead, is passed to the analog electronic domain in the form of back e.m.f. voltage amplitude according to: $V_{backe.m.f.}=k_{e.m.f.}dw$, where $k_{e.m.f.}$ represents a numerically estimated constant: $k_{e.m.f.}=1.02\times {10}^{-12\frac{\mathrm{Vs}}{\mathrm{mm}}}$.

Figure 8 shows an implemented numerical model of the equivalent MEMS device developed through ESRF Radia^® package and solved through the 3D Finite Volume Method (3DFVM). Since most of the constitutive parts do not have any ferromagnetic properties and, consequently, do not interact during the simulation, the system was modeled through a reduced number of components: the planar micro-coil, the permanent magnet, and a fictitious element depicting the PDMS diaphragm. The coil was grounded, while the fictitious disk together with the magnet were continuously updated in position to have a correspondence between the magnetic model and the MBD model.

Due to the lack of iron-based materials, the part that required the most computational effort and, consequently, required a volume partition was the NdFeB magnet. The best compromise in solution convergence and computing performances was obtained with 4 × 4 × 4 volume subdivisions along the X-axis, Y-axis, and Z-axis, respectively. To best fit the cylindrical geometries of the components, the polyhedron finite volume shapes were adopted. Once the 3D volume objects were created, the materials were applied and the interaction matrix was obtained, the relaxation procedure was executed and all the magnetization vectors were extracted by the single tiny volumes, allowing us to obtain the magnetic field intensity at each point of the 3D space.

The GUI circuit schematic of the Spice-based equivalent model is reported in Figure 9. The circuit was implemented for the NGSpice engine simulator. The planar micro-coil was modeled here through a parallel between the stray capacitance and the series of the AC resistance, the equivalent inductance, and the back e.m.f. generator; the constitutive parameters values were updated for every time step call of the simulation based on the solved variables passed back by the leading Simulink routine. The circuit was fed by a $10\mathsf{\Omega }$ source generator and a switch. In order to replicate a real case application and to enhance the potentiality of the proposed architecture, an active actuation control was introduced; this microcontroller replica was implemented through an equivalent generator and a MOSFET. The circulating current was measured through the voltage drop across the shunt resistor terminals.

Finally, the mechanical model of the MEMS actuator was developed and simulated as a multi-body (MBD) model. The equivalent MBD model was constituted by three main sub-assemblies:

• The NdFeB permanent magnet (red body);
• The PDMS diaphragm (blue body);
• The frame (green body) was composed by the PMMA plate, the micro-coil, the polyimide, and the glass substrate.

In the following are described the modeling techniques chosen for the replication of the kinematic and dynamic properties characteristics of an MEMS electromagnetic actuator and are schematized in Figure 10.

• The kinematic deformation properties of the membrane were modeled through a prismatic joint ➀ between the fictitious diaphragm disk and the frame.
• The magnet was kinematically constrained through a fixed joint ➁ to the fictitious membrane.
• A fixed constraint ➂ settled the frame to the ground.
• A magnetic force ➃ was applied between the magnet and the frame. Its absolute value and the verse of application represented one of the co-simulation state variables exchanged among the different simulation platforms.
• A spring-damper equivalent force ➄ was added to replicate the elastic and dissipative properties of the PDMS diaphragm. The constitutive parameters values were estimated according to Equations (12)–(15).
• The measure $w$ and $dw$ of the equivalent deflection was measured between the frame and the disk and was passed back to the leading routine as presented in Figure 6 to compute the induced e.m.f.

Obviously, from the point of view of multibody simulation, the reported model represents a very simplified application. The case in analysis is to be understood as a basic explanatory model for the complete description and validation of the presented multidomain co-simulation architecture.

5.5. Results Comparison and Validation

The validation of the proposed co-simulation architecture was verified by comparing the obtained numerical results with the experimental measurements reported in the referred article [15].

The requested simulation time for the single analysis was around 35 min with an 11th Gen Intel Core i7 CPU, NVIDIA GeForce GTX 1650 GPU and DDR4 32 Gb RAM.

Figure 11 depicts on two different graphs (a) the laser-based displacement measure obtained at different levels of applied current and (b) the numerical trends of the estimated deflection amplitude in a 0.8 s time window for a circulating current equal to: [0.2; 1.0] A. For each one of the reported configurations, the trends show an initial increasing monotone transition phase and, subsequently, a maintenance phase into a deformed configuration. The final maintenance values were very similar to what was reported by the experimental measurements. It is important to underline that the greatest deviation between measurements and numerical estimations occurred for current values equal to $1.0\mathrm{A}$. This response could come from non-linear elastic behavior obtained for such large strain deformation.

To emphasize the capabilities and the wide range of simulation possibilities provided by the proposed multi-domain architecture, a series of results obtained from the previously described actuator model are reported below. In detail, we studied the dynamic response of the device through a variable frequency-modulated excitation. The signal provided by the microcontroller equivalent generator (see Figure 12) commanded the MOSFET gate according to a PWM voltage characterized by three main frequencies ([80; 90; 100] Hz) and duty cycles ([50%; 60%; 70%]).

Figure 13 shows the dynamic response of the device. It is important to notice how the estimated deflection strictly follows the PWM frequency changes occurring at intervals of 200 ms. After 600 ms, the excitation signal ended, and the device returned to the quiescent configuration through a damped oscillatory motion.

Finally, Figure 14 and Figure 15, respectively, show some of the most interesting results obtained from the magnetic and electronic domains.

In detail, Figure 14 depicts the characteristic trend of the electro-magnetic actuation force. The first graph shows the trend during the entire actuation window, while the lower one shows a detailed magnification of the characteristic waveform. Figure 15, instead, represents the time history of the voltage measured on the MOSFET drain in a time interval of [0.18; 0.22] s, in which are clearly visible the under- and over-shoot (see the graph below) characteristics of a real MOSFET Spice model.

The reported series of results highlights how, although on a simple but extremely representative case study, the presented co-simulation methodology allows researchers to investigate in depth the strongly coupled phenomena that characterize general electro-magnet-mechanical systems actuated by a firmware control logic.

6. Conclusions

In this study, a novel time-based multi-domain architecture is presented, and its capabilities and accuracy are verified. The method uses a multi-physics approach to embed under the same simulation platform general purpose electro-magneto-mechanical problems. Two novel Python-Matlab based bi-directional communication routines are presented, with the aim of integrating two consolidated simulation tools and languages for the electronic and magnetic domains in multi-physics models: a Spice-based analysis and a 3D Radia-based finite volume methodology, respectively. It was established that for the magnetostatic analysis the FVM principle presents significant improvements both in reliability and performance compared with the classic FEM analysis. The integration of the Spice-based electronics aims to extend this design boost to those circuits, which, due to the application complexity, require a direct interaction between different physical domains, such as in the case of mechatronic applications.

It is important to underline how the proposed co-simulation algorithm allows the interaction of the different phenomena, while keeping the characteristic frequencies of the studied events and the properties of the called integrators fully decoupled, so as to optimize the computational load of the CPU without incurring aliasing phenomena.

Finally, a detailed validation is faced through a complete comparison of quasi-static measurements [15] obtained from an MEMS electromagnetic actuator. The dynamic regime analysis is also considered. To enhance the capabilities and the reliability of the proposed architecture, a real case scenario is studied. The system modeled with all the equivalent constitutive components is connected to a fictitious microcontroller equipped with a firmware control logic which, through different PWM signals, allows us to study the full dynamic response. In detail, signals with exciting frequencies of [80; 90; 100] Hz and duty cycles equal to [50%; 60%; 70%] are analyzed.