Some Case Studies in Environmental and Industrial Mathematics

Abstract

This presentation deals with four case studies in environmental and industrial mathematics developed by the mathematical engineering research group (mat+i) from the University of Santiago de Compostela and the Technological Institute for Industrial Mathematics (ITMATI). The first case involves environmental fluid mechanics: optimizing the location of submarine outfalls on the coast. This work, related to shallow water equations with variable depth, led us to develop a theory for numerical treatment of source terms in nonlinear first order hyperbolic balance laws. More recently, these techniques have been applied to solve Euler equations with source terms arising from numerical simulation of gas transportation networks when topography via gravity force is considered in the model. The last two problems concerns electromagnetism. One of them is related to nondestructive testing of car parts by using magnetic nanoparticles (the so-called magnetic particle inspection, MPI): mathematical modelling of magnetic hysteresis to simulate demagnetization. Finally, we present a mathematical procedure to reduce the computing time needed to achieve the stationary state of an induction electric machine when using transient numerical simulation.

1 Introduction

Four case studies developed by the Research Group in Mathematical Engineering from the University of Santiago de Compostela (USC) and the Technological Institute for Industrial Mathematics (ITMATI) are considered. Two of them are related to fluid mechanics. The first one was developed in the framework of a contract with the Ministry of Public Works of Galicia and concerns shallow water flows in a domain with variable depth. The second one deals with gas flow in transport networks and has been done for the Reganosa company. From the mathematical point of view both are modelled with systems of nonlinear hyperbolic partial differential equations with source terms and the goal is to set up suitable finite volume discretization of the source terms.

The other two case studies concern electromagnetism. The goal of the first one, that has been financed by CIE Automotive company, is numerical simulation of magnetization and demagnetization processes in magnetic particle inspection procedures. Finally, the last case study is related to numerical solution of electric machines with optimal design in view. The underlying mathematical problems are, respectively, mathematical and numerical analysis of models for electromagnetic hysteresis, and methods to determine appropriate initial conditions for transient electromagnetic simulations, in order to attain the steady state as soon as possible.

2 Environmental Flows. The Shallow Water Equations

The technical goal of this work, commissioned by the Galician government to our research team in the eighties, was to determine the optimal location of submarine outfalls along the coast of the Galician rias. For this purpose several steps were done involving modelling, simulation and optimal control:

• To compute the velocity field due to tidal currents and wind which was done by using the shallow water equations

• To solve a mathematical model giving the evolution and dispersion of some pollution indicators as fecal coliforms or biochemical oxygen demand (BOD)

• To formulate and solve some constrained optimal control problems related to outfall position and management of wastewater treatment systems.

Regarding the first step, as the shallow water equation is a nonlinear system of hyperbolic partial differential equations, numerical methods developed in the eighties of the last century for Euler equations can be applied to its numerical solution. We mean finite volume methods combined with approximate Riemann solvers. The unexpected problem we found was related to the discretization of the source term which is present in the shallow water equations when the bottom is not flat. In order to give some insight we refer to Fig. 1: we have solved the shallow water equations by using a finite volume scheme with the Van Leer Q-scheme as approximate Riemann solver for flux term upwinding, and a centred scheme to discretize the source term arising from non-flat bottom. We have considered a static configuration in a closed channel, more precisely, the initial condition (and then the solution along the time) corresponds to water at rest. In the left plot one can see the computed water level which is a quite good approximation. However, the right plot shows the computed velocity which varies between around $-60$ and 80 m/s while the exact velocity is null.

Fig. 1

Shallow water. Centred discretization of the source term. Computed water level (left) and computed velocity (right). Notice that the zero line is the result of a numerical simulation using [10]

Motivated by this problem, in the old paper [10] we developed a general methodology to discretize source terms in nonlinear systems of first order hyperbolic partial differential equations. In particular, our methods solve exactly the previous static problem. This paper is considered a seminal work in the theory of well-balanced schemes for numerical solution of conservation laws with source terms, an active field of research during the last years. Moreover, thirty years later, this methodology was applied by our research group to a different problem: Euler equations with gravity, more specifically, to numerical simulation of gas transportation networks on non-flat topography.

3 Gas Network Simulation

This industrial demand from the Reganosa company consisted in writing a software code for transient numerical simulation of a gas transport network. In Fig. 2 the high-pressure Spanish gas network is shown. Besides the great number of pipes, it includes entry (emission) and exit (consumption) points, underground storages and, more importantly, compression stations. The latter are needed to compensate the pressure drop along the network due to viscous friction of the gas on the pipe walls.

Fig. 2

Spanish gas transport network

3.1 Mathematical Modelling: Homogeneous Gas Flow in a Pipe

The mathematical model for gas flow in a pipe consists of Navier-Stokes equations for compressible flows. More precisely, it involves the mass, momentum and energy conservation laws and some additional equations: the state equations for real gases and the Darcy-Weisbach law for turbulent friction between gas and pipe walls combined with Colebrook equation to compute the friction factor. As the pipe length is much larger than the area of its cross-section we can use a 1D model:

$$\begin{aligned} \frac{\partial {\rho }}{\partial t}(x,t) +\frac{\partial ({\rho } {v})}{\partial x}(x,t)=&0,\\ \frac{\partial ({\rho } {v})}{\partial t}(x,t)+\frac{\partial ({\rho } {v}^2+p)}{\partial x}(x,t)=&\underbrace{-\frac{\lambda {\rho }(x,t)}{2D}|{v}(x,t)|{v}(x,t)}_{\text{ friction }} \underbrace{ -g{\rho }(x,t)h^{\prime }(x),}_{\text{ gravity } \text{ force }}\\ \frac{\partial ({\rho } {E})}{\partial t}(x,t)+\frac{\partial (({\rho }{E}+p) {v})}{\partial x}(x,t)=&\underbrace{{-g {\rho }(x,t) v(x,t) h^{\prime }(x)}}_{\text{ power } \text{ of } \text{ gravity } \text{ force }}\\&+\ \underbrace{{\alpha \frac{4}{D} (\theta _{ext}(x,t)-\theta (x,t)).}}_{\text{ heat } \text{ exchange }}\\ \text{ Thermodynamic } \text{ equation } \text{ of } \text{ state: }\, p=&Z(\theta ,p)\rho R\theta \\ \text{ Caloric } \text{ equation } \text{ of } \text{ state: }\, e=&E-\frac{1}{2}|v|^2 \quad \text{ with } \\ e=&\ \hat{e}(\theta )=\hat{e}(\theta _0)+\int \limits _{\theta _0}^\theta c_v(s)\,ds \end{aligned}$$

• $\theta$ is absolute temperature (K)

• p is pressure (Pa)

• $Z(\theta , p)$ is the compressibility factor (dimensionless)

• E is the specific total energy (J/kg)

• e is the specific internal energy (J/kg)

• $\theta _0$ is a reference temperature (K)

• $c_v(\theta )$ is the specific heat at constant volume (J/(kg K)).

3.2 Numerical Solution: One Single Pipe with Homogeneous Gas

Numerical methodology for solving the compressible Euler equations for homogeneous mixtures of perfect gases without sources has been well established since the eighties of the last century. For instance, one can use a simple first-order method consisting of Euler explicit for time discretization, finite volume method for space discretization, and approximate Riemann solvers (e.g., van Leer’s Q-Scheme) for upwind discretization of the flux term (see, for instance, [24]). However, when source terms are present (e.g., the gravity term with variable heigth), numerics is more difficult and similar to the shallow water equations the use of well-balanced schemes is mandatory. This means that the discretization of source terms also needs some upwinding. In the last years many papers devoted to numerical solution of Euler equations with gravity have been written. Let us mention, for instance, [13,14,15, 23, 25, 27].

In order to highlight the need of using an upwind discretization of the source terms, we consider the following very simple test problem: h(x) in the gravity source term is an arbitrary function and we look for a static isothermal solution, i.e., satisfying $v(x)=0,\quad \theta (x)=\theta _{ext},\quad \forall x\in (0,L).$ It is easy to see that the exact solution is given by $v(x)=0$, $\rho (x) =\rho _0 \exp \left( \displaystyle -\frac{g}{R \theta _{ext} } \big (h(x)-h_0\big )\right) ,$ and $p(x) =R \theta _{ext} \rho _0 \exp \big (-\displaystyle \frac{g}{R \theta _{ext} } \big (h(x)-h_0\big )\big ).$ For the data given in Table 1, the computed mass flow rate is shown in Fig. 3 as well as the exact solution which is null. One can see that the former is very bad, oscillating between around −10 and 10.

Table 1 Data for static isothermal test

Fig. 3

Mass flux, (kg/(m$^2$ s)). Computed with centred discretization of source terms (black) and exact (red). The horizontal axis is the distance to the origin of the pipe

By using the general methodology developed in [10], we have proposed a discretization of the gravity term in [7] leading to a well-balanced scheme that reproduces the null solution exactly.

3.3 Network with Heterogeneous Gas

Simulation of heterogeneous gas flowing in a network is more difficult. New problems arise: junction modelling, gas quality simulation. These issues have been addressed in papers [8, 9].

Fig. 4

The Reganosa network (Galicia. Spain)

3.4 Experimental Validation in a Real Small Network

The code has been used for a small gas network and the results have been compared to real measurements. The network can be seen in Fig. 4.

Topography is quite irregular as can be seen in Fig. 5. Results and measurements corresponding to mass flow rate and pressure for some particular nodes are shown in Figs. 6 and 7, respectively.

Fig. 5

Height function for edge #9

Fig. 6

Mass flow rate at node 01A. Blue: real measurement. Red: computed with a homogeneous gas model. Green: computed with a heterogeneous gas model

Fig. 7

Pressure at node I-015. Blue: real measurement. Red: computed with a homogeneous gas model. Green: computed with a heterogeneous gas model

4 Non-destructive Testing: Magnetic Particle Inspection (MPI)

MPI is a non-destructive testing technique to detect near-surface defects in ferromagnetic pieces. The process is as follows: firstly, the workpiece is magnetized. Then, the presence of a surface discontinuity in the material allows the magnetic flux to leak, since air cannot support as much magnetic field per unit volume as metals. In order to identify a leak, ferrous particles, either dry or in a wet suspension, are applied to the workpiece. Then they are attracted to an area of flux leakage and form what is called an indication which is evaluated to determine its nature. Since cracks are more easily detected when they are perpendicular to the induced field, two magnetizations are made: circular and longitudinal. After inspection, a final demagnetization step is required for subsequent processing of the workpiece. In the next subsection we introduce an axisymmetric model for circular magnetization and present some numerical results (Figs. 8 and 9). Further details can be found in Refs. [2, 4,5,6].

Fig. 8

Magnetic particle inspection

Fig. 9

Crack indication. Circular magnetization. Longitudinal magnetization

Fig. 10

Circular magnetization

4.1 Circular Magnetization. Axisymmetric Model

Let us introduce a mathematical model for circular magnetization. Thanks to axisymmetry, it can be written on a meridional section (see Fig. 10).

Given I(t), the magnetizing or demagnetizing current, and an initial condition $H_0$, find $H_{\theta }$ in $\varOmega \times (t_0,T]$ such that

$$\begin{aligned} \frac{\partial B_\theta }{\partial t}\boldsymbol{e}_{\theta } + \mathop {\mathbf {curl}}\left( \frac{1}{\sigma } \mathop {\mathbf {curl}} (H_{\theta }\boldsymbol{e}_{\theta })\right) = \boldsymbol{0}&\quad \text{ in } \varOmega \times (t_0,T],\\ H_{\theta }(0,z,t) = 0&\quad \text{ on } (0,L)\times (t_0,T],\\ H_{\theta }(R_{\mathcal {S}}(z),z,t)= \frac{I(t)}{2\pi R_{\mathcal {S}}(z)}&\quad \text{ on } (0,L)\times (t_0,T],\\ \frac{\partial H_{\theta }}{\partial z}(\rho ,z,t) = 0&\quad \text{ on } \left( \varGamma _1\cup \varGamma _2\right) \times (t_0,T],\\ H_{\theta }(\rho ,z,t_0) =H_0(\rho ,z)&\quad \text{ in } \varOmega . \end{aligned}$$

and

$$ B_\theta (x,t)=\mathcal {B}(H_{\theta }(x,.),\xi (x))(t),$$

where $\mathcal {B}$ is a scalar hysteresis operator to be defined later.

4.2 Hysteresis Modelling

Mathematical modelling of hysteresis is now a well established subject (see, for instance, the reference books [11, 12, 17,18,19, 26]). Let us summarize the main issues of the theory. We consider a system whose state is characterized by two scalar variables, u and w, which are assumed to depend continuously on time t. In our case $u=H_\theta$ and $w=B_\theta$. The value of w(t) is determined by u(t) and by the values of $u(\tau )$ for $\tau < t$. Let us introduce some basic definitions and notations (Fig. 11).

Fig. 11

Hysteresis major and minor loops

At any instant t, w(t) depends on the previous evolution of u, and on an initial state of the system to be called $\xi$. We can formalize this as follows:

$$ w(t) = \mathcal {F}(u,\xi )(t) \quad \forall \, t \in [0, T]. $$

Here $\mathcal {F}(\cdot ,\xi )$ represents an operator between suitable spaces of time-dependent functions. Notice that $\mathcal {F}$ is non-local in time. A particular example of hysteresis operator is the Preisach operator:

$$\begin{aligned}&\mathcal {F}: C^0([0,T])\times Y \longrightarrow C^0([0,T]), \\&[\mathcal {F}(u,\xi )](t):=\int \limits _{\mathcal {T}}[h_\rho (u,\xi (\rho ))](t)p(\rho )d\rho , \end{aligned}$$

where $\mathcal {T}$ is the Preisach triangle, $0<p\in L^1(\mathcal {T})$ is the Preisach function which is determined by physical experiments for each material (see Fig. 12), $h_\rho$ is the relay function (see Fig. 13) and $\xi :\mathcal {T} \rightarrow \{-1, 1\}$ is a Borel measure representing the initial magnetic state.

Fig. 12

Preisach triangle (left) and an example of Preisach function (right)

The classical Preisach model is built with the so-called rate-independent relay: let us fix any pair $\rho := (\rho _1, \rho _2) \in \mathbb {R}^2,\, \rho _1 < \rho _2$. For any continuous function $u : [0,T] \rightarrow \mathbb {R}$ and any $\xi \in \{-1, 1\}$, we define $h_\rho (u,\xi )$ as follows.

Let $t_1<\ldots<t_N<t$ be such that $u(t_i)\in \{\rho _1,\rho _2\}$. If $\{t_i\}=\emptyset$ or $t=0$, then

$$ h_\rho (u,\xi )(t)=\left\{ \begin{array}{cl} -1 &{} \mathrm {if }\, u(t)\!\le \rho _1, \\ \xi \!\!\! &{} \mathrm {if }\, \rho _1\!<\!u(t)\!< \rho _2, \\ 1 \!\!\!&{} \mathrm {if}\, u(t)\!\ge \rho _2, \\ \end{array}\right. $$

else

$$ h_\rho (u,\xi )(t):=\left\{ \begin{array}{cl} 1 &{} \mathrm {if}\,\, u(t_N)=\rho _2, \\ -1 &{} \mathrm {if}\,\, u(t_N)=\rho _1. \end{array}\right. $$

If we split $\mathcal {T}=S_u^{+}(t)\cup S_u^{-}(t)$, where

$$ S_u^{\pm }(t)=\left\{ (\rho _1,\rho _2)\in \mathcal {T}: [r_{\rho }(u,\xi )](t)=\pm 1\right\} , $$

then

$$ [\mathcal {F}(u,\xi )](t):=\int \limits _{S_u^{+}(t)}p(\rho )d\rho -\int \limits _{S_u^{-}(t)}p(\rho )d\rho . $$

Fig. 13

Classical relay operator

Fig. 14

Input u(t) (left) and its corresponding splitting of Preisach triangle (right)

We present some results obtained by solving the above model for a real crankshaft (see Fig. 14 for input data). Figure 15 shows the remanent magnetization after the circular magnetization process. In its turn, Fig. 16 shows the applied demagnetization current and the remanent magnetization after demagnetizing.

Fig. 15

Remanent magnetization

Fig. 16

Demagnetization current (left) and remanent magnetization after demagnetizing (right)

5 Accelerated Simulation of Electric Machines

In the design of electric machines (see Fig. 17), numerical simulation is an important tool. The engineer needs to know the behaviour of the machine in steady regime. In particular, he/she wants to know the torque. In order to get this steady state, finite element methods are used to solve a transient nonlinear system of PDEs derived from Maxwell equations, coupled with electrical circuit equations, starting from an (arbitrary) initial condition until the steady state is achieved. The time for this transient model to attain the steady state highly depends on the choice of the initial condition. When an unappropriate value is prescribed (for instance, when it is set to zero), a very long CPU time is needed to reach the steady state solution. Therefore, techniques leading to a suitable initial condition are in high demand and in the literature we can find several approaches to the problem. Let us mention, for instance, time periodic finite element methods [21], time periodic-explicit error correction methods [16], time differential correction [20], parareal algorithms [22]. A common drawback for these methods is the need of choosing a suitable time interval in which the solution is assumed to be periodic: the so-called effective period. Indeed, magnetic fields in rotor and stator oscillate at different frequencies and the common time at which both are periodic is generally quite large. However, the periodicity condition has to be defined in a short time interval for the method to be useful. Our methodology aims to compute a suitable initial condition and has the advantage of making use of periodicity property only in the rotor bars, so the above limitation does not apply. Moreover, the computational cost of our approach does not depend on the size of this period, and the number of unknowns is very small in comparison with the previously mentioned methods.

This work has been developed under contract with the company Robert Bosch GmbH from Stuggart (Stefan Kurz, Marcus Alexander). It has given rise to a Spanish patent. A detailed description of the methodology has been published in papers [1] and [3].

Fig. 17

From Wikimedia Commons by Mtodorov 69 under license CC-BY-SA-3.0

Main parts integrating an induction motor.

Fig. 18

A quarter of the geometric domain at time $t=0$ (left) and $t>0$ (right). Modification of a picture provided by Robert Bosch GmbH

5.1 Description of the New Methodology

The main lines of the developed methodology can be described for a toy model. Let us consider a simple series circuit with an inductor and a resistor,

$$\begin{aligned} \mathrm {L}\dot{\mathrm {I}}(t) + \mathrm R\mathrm {I}(t) =\mathrm E(t), \end{aligned}$$

with the electromotive force

$$ \mathrm E(t)=\mathbb {E}\sin (\omega t) $$

The general solution is

$$\begin{aligned} \mathrm {I}(t)=\underbrace{A\mathrm e^{-\frac{\mathrm R}{\mathrm {L}}t}}_{\text{ transient } \text{ part }}+\underbrace{\frac{\mathbb {E}}{|\mathcal {Z}(\omega )|}\sin (\omega t-\varphi (\omega ))}_{\text{ steady } \text{ solution }} \end{aligned}$$

where $\mathcal {Z}(\omega )=\mathrm R+\omega \mathrm {L}i\in \mathbb {C}$ is the impedance of the circuit and $\varphi (\omega )$ its argument. We have two opposite extreme situations:

• If $\frac{\mathrm RT}{\mathrm {L}}\gg 1$, then the exponential vanishes quickly independently of the initial condition

• If $\frac{\mathrm RT}{\mathrm {L}}\ll 1$ then the transient part strongly depends on the initial condition. Moreover, in this case

  $$ \varphi (\omega ) \approx \frac{\pi }{2} \text{ and } |\mathcal {Z}(\omega )|\approx \omega \mathrm {L}$$

  and hence

  $$ \mathrm {I}(t) \approx A \mathrm e^{-\frac{\mathrm R}{\mathrm {L}} t} +\frac{\mathbb {E}}{\omega \mathrm {L}} \cos (\omega t ). $$

If the equation is solved for $\mathrm {I}(0) = 0$, then the solution is approximately given by

$$ \mathrm {I}(t) \approx -\frac{\mathbb {E}}{\omega \mathrm {L}} \mathrm e^{-\frac{\mathrm R}{\mathrm {L}} t}+ \frac{\mathbb {E}}{\omega \mathrm {L}} \cos \omega t, $$

so it includes a transient part. However, if the equation is solved for

$$ \mathrm {I}(0) = \frac{\mathbb {E}}{\omega \mathrm {L}}. $$

then $A=0$ and the transient part is close to zero from the beginning. The important remark is that, if $\frac{\mathrm RT}{\mathrm {L}}\ll 1$ then the above initial condition can be obtained without solving the ODE, as follows:

• Firstly, the term involving the resistor can be neglected

• Then, we integrate the equation twice: first between 0 and t and then between 0 and T. We get

  $$ \mathrm {L}\int \limits _0^T\mathrm {I}(t)\,dt-\mathrm {L}T\mathrm {I}(0)=\mathbb {E}\int \limits _0^T(T-s)\sin \omega s\,ds=\frac{\mathbb {E}T}{\omega } $$

• Moreover, since the steady solution is harmonic then $\int \nolimits _0^T\mathrm {I}(t)\,dt=0$ and from the above equation we deduce

  $$ \mathrm {I}(0)=\frac{1}{\mathrm {L}T}\frac{\mathbb {E}T}{\omega }=\frac{\mathbb {E}}{\omega \mathrm {L}} $$

which is the suitable initial condition previously obtained. The interesting feature of this method is that it can be used in more general settings; in particular, to the model of induction machines with squirrel cage. In this case, the problem to be solved is the following:

Given currents along the coil sides $I_n(t)$, $n=N_b+1,\ldots ,N_c$, and initial currents along the bars $y^0_n$, $n=1,\ldots ,N_b$, find, for every $t\in [0,T]$, currents $y_n(t)$, $n=1,\ldots ,N_b$, along the bars such that $y_n(0)=y^0_n$, $n=1,\ldots ,N_b$, and

$$\begin{aligned} {\mathcal R}^b\frac{d}{dt} \mathbf {\mathcal {F}}\left( t,\mathbf {y}^{\,b}(t)\right) +\left( {\mathcal R}^b+\left( {\mathcal A}^b\right) ^{\mathrm {T}}{\mathcal B}^{-1}\left( {\mathcal A}^b\right) \right) \mathbf {y}^{\,b}(t)+ \lambda (t)\left( {\mathcal A}^b\right) ^{\mathrm {T}}\left( \begin{array}{c}\mathbf {0}\\ \mathbf {e}\end{array}\right)&=\mathbf {0},\\ {\mathcal A}^b\mathbf {y}^{\,b}(t) \cdot \left( \begin{array}{c}\mathbf {0}\\ \mathbf {e}\end{array}\right)&= 0, \end{aligned}$$

where $\mathbf {\mathcal {F}}:[0,T]\times \mathbb {R}^{N_b}\longrightarrow \mathbb {R}^{N_b}$ is the nonlinear operator defined as

$$\begin{aligned} \mathbf {\mathcal {F}}(t,\mathbf {w}):=\left( \int \limits _{\varOmega _1} \sigma A(x,y,t)\ dx \, dy, \ldots , \int \limits _{\varOmega _{N_b}} \sigma A(x,y,t)\ dx \, dy \right) ^{\mathrm {T}}\in \mathbb {R}^{N_b}, \end{aligned}$$

for $\quad t\in [0,T],\mathbf {w}\in \mathbb {R}^{N_b},$ with A(x, y, t) the solution to the following nonlinear magnetostatic problem:

Given a fixed $t\in [0,T]$, currents along the coil sides $I_n(t)$, $n=N_b+1,\ldots ,N_c$, and $\mathbf {w}\in \mathbb {R}^{N_b}$, find a field A(x, y, t) such that

$$\begin{aligned} - \mathop {\mathrm {div}}\nolimits (\nu _0\, \mathbf {grad}A)&= 0 \quad \text{ in } \varOmega _{0}^{\mathrm {rot}}\cup r_t\left( \varOmega _{0}^{\mathrm {sta}}\right) ,\\ - \mathop {\mathrm {div}}\nolimits (\nu _0\, \mathbf {grad}A)&=\frac{w_n}{\mathop {\mathrm {meas}}\nolimits (\varOmega _n)} \quad \text{ in } \varOmega _n, \, n=1,\ldots ,N_b,\\ - \mathop {\mathrm {div}}\nolimits (\nu _0\, \mathbf {grad}A)&=\frac{I_n(t)}{\mathop {\mathrm {meas}}\nolimits (\varOmega _n)} \quad \text{ in } r_t(\varOmega _n), \, n=N_b+1,\ldots ,N_c,\\ - \mathop {\mathrm {div}}\nolimits (\nu (\cdot ,|\mathbf {grad}A|)\, \mathbf {grad}A)&=0 \quad \text{ in } \varOmega _{\mathrm {nl}}^{\mathrm {rot}}\cup r_t\left( \varOmega _{\mathrm {nl}}^{\mathrm {sta}}\right) , \end{aligned}$$

with suitable transmission and boundary conditions.

5.2 Numerical Experiments with Real Electric Machines

We present the numerical results obtained for a particular induction machine with squirrel cage rotor. Firstly, we use our method to get a suitable initial condition. Next, we solve the transient model with this initial condition and compare the time needed to reach the steady-state with the one needed by taking null initial condition. The electric machine we have used for numerical experiments can be seen in Figs. 18 and 19. For confidentially issues it is a modification of a picture provided by Robert Bosch GmbH. Red, yellow and blue colors correspond to the three different phases. It is composed by 36 slots in the rotor and 48 slots in the stator. It is a three-phase machine having 2 pole pairs with 12 slots per pole. The source currents are characterized by an electrical frequency $f_c$ and a RMS current $I_c$ through each slot. The currents corresponding to each phase of the stator are defined as

$$\begin{aligned} I_A(t)&= \sqrt{2} \,I_c \cos \left( 2 \pi f_c t\right) , \\ I_B(t)&= \sqrt{2} \,I_c \cos \left( 2 \pi f_c t + \frac{2 \pi }{3}\right) , \\ I_C(t)&= \sqrt{2} \,I_c \cos \left( 2 \pi f_c t - \frac{2 \pi }{3}\right) . \end{aligned}$$

We have considered four operating points corresponding to different electrical sources in the stator and different rotor velocities. They are described in Table 2. The physical time to reach the steady state for the different operating points can be seen in Table 3. Finally, in Fig. 20, the computed torque and current along the transient simulation are shown for operation point # 4.

Fig. 19

Transversal section of an induction electric motor with squirrel cage

Table 2 Operation points for numerical tests

Table 3 Time to get the steady state with null initial condition and with the one obtained by the new method

Fig. 20

Op. Point 4. Torque versus time (left). Current in bar 1 versus time (right)

Notes and Comments.

• We have presented four case studies in industrial mathematics, all related with numerical simulation by partial differential equations

• In addition to the industrial outcome, in all cases scientific papers related to the developed methods have been published

• This shows that industrial problems usually lead to new mathematical developments

• Industrial mathematics is a nice area with good opportunities for young mathematicians willing also to learn other scientific disciplines

• Postgraduate studies mixing applied mathematics and areas of application as physics, chemistry, biology, medicine, economy, etc. are a good initial step to develop a career in this promising area of increasing interest for companies and research institutions.