1 The Intention of the Study

The main objective of the paper is to develop a program to solve electromagnetic problems with numerical solutions using the link between mathematical/geometrical operations and physical phenomena. The importance of writing and developing a program is emphasized to understand the structure of these complex equations, learn the properties of the solution method, and conduct research on the method. The study is realized on a specific ferromagnetic domain with different scenarios for educational purposes.

2 Introduction

The solution of real magnetic circuits using the 2D analytical solution is time-consuming and very complex for electrical engineering, researchers, students, and even for today's powerful computers because most magnetic circuits are solved with material nonlinearity. Therefore, numerical methods, such as the finite element method (FEM) or the finite difference method (FDM), which transform a continuous domain into a finite partition, are used to create a simplified and understandable model to solve magnetic circuits or Maxwell equations [1, 2]. Since numerical modeling in electromagnetism is a tool to solve rigor and complex problems, the use of software packages for modeling and simulation is extremely necessary. Magnetic field analyses have been studied for various problem types with different methodologies, which deal with electromagnetic field techniques in detail with their strengths and weaknesses. To overcome complex and nonlinear numerical equations in electromagnetic problems, some studies have preferred simulation programs with simplified models using 2D axisymmetry [3,4,5]. An electromagnetic actuator is studied to analyze heating electromagnetic field distribution and decrease the power losses using the FEMM simulation program [6]. In [7], some discretization techniques based on two (or more) FEM subspaces for solving partial differential equations were presented. It proposed a new technique to solve a large class of nonlinear elliptic boundary value problems. FEM and FDM have been compared in [8, 9], and the advantages of FDM over FEM in some cases are discussed. Electromagnetic numerical problems indeed have complex, strongly nonlinear, and multimodal solution algorithms, and thus, the lecturer motivates the students to recognize that basic electromagnetic knowledge is an open field for research using powerful mathematical software and coding. The software or coding is characterized by the proximity between the way to define discrete problems and their symbolic mathematical expressions written in a mathematical form that is adapted to a chosen numerical method, directly constituting data through structured computational objects. This allows the theory and mathematical bases of numerical methods to be directly followed by their practical applications [10, 11]. The vector computation of electromagnetic fields presents a mathematical stumbling set to many students, which discourages their understanding of physical concepts. The present courses are mainly based on the behavior of electromagnetic fields upon differential and integral vector calculus. Because the visualization of vector fields in the students' minds is challenging, the fundamentals of electromagnetics courses often become exercises where students only memorize a set of formulas, making minimal connections between the physical phenomena and mathematics. If the formulation of electromagnetic laws in numeric form should be given to beginning students directly with computer-based solutions, the attractiveness of the subject can be achieved in real-time, and scalar calculations displayed with color graphics can provide a powerful learning environment [12, 13]. However, except for a few education coding programs, most research is now conducted using sophisticated simulation programs.

In this study, Maxwell's equations, which are the basic equations of electromagnetic fields, are solved with simplified 2D static-elliptical electromagnetic problem analysis coding and software tools for the numerical solution of a physical problem. To accomplish this, the problem domain is first created and visualized, and then the vector and scalar potential formulations are applied to the specified domain with boundary conditions. Finally, the finite element magnetic method (FEMM) simulation is investigated in the same region to give an analogy and compare the developed model results. The coding for building region and related formulations are given in the text and Appendix. The coding method meets an educational need by incorporating the complex Maxwell equations platform for teaching electromagnetic engineering. Vector potential formulation can combine all the conditions that must be satisfied into a single equation. If the vector potential is discovered in this type of equation, differentiating \(A\)(vector potential) yields \({\text{\rm B}}\)(magnetic field density) and \(H\)(magnetic field intensity). Many tools and software have been developed over the years to solve this problem [9,10,11]. Mathematica software was chosen for coding because of its powerful symbolic language, intelligent automation from algorithm selection to plot layout, programming capabilities, and intuitive user interface design [12, 15]. The approach used is suitable for an undergraduate and a graduate student, or even for a senior engineer, depending on the level of depth and rigor.

3 Computational Tools to Solve Electromagnetic Problem

3.1 Physical Structure of the Magnetic Device

The physical structure, developed in Mathematica, of the problem in dimensions being studied is depicted in Fig. 1. The magnetomotive force (mmf) is an ampere multiplied by the number of turns in the coil \((N\,x\,\,I)\); hence, the required mmf can be easily determined. The ferromagnetic material properties and the coil current are the main characteristics used to determine the magnetic field of a ferromagnetic device. The numerical solution of the Maxwell equation is performed using the mmf, the magnetic structure dimensions, and the boundary conditions. There are two permeability effects shown in blue and red in Fig. 1. Additionally, they are analyzed in vector and scalar potential formulations using Eqs. (1) and (2). The appropriate boundary conditions for the structure are determined using Eqs. (3) and (4). The results are evaluated by matching Eqs. (1) and (2) with Eqs. (3) and (4) in different ways, and the results of the matching options remained consistent in all cases. The problem solution is realized in two regions: one is the entire region, including the gray part in Fig. 1 known as the “domain,” and the other is the only mmf-produced region, which is depicted by blue and red parts together and is dubbed the “solid.” Another region, known as the “coil” region, produces mmf and is colored black.

Fig. 1
figure 1

Schematic diagram of the magnetic structure (solid), dimensions are in centimeters

Full size image

The parameter values of the studied ferromagnetic structure with different materials are shown in Table 1. In the table, the relative permeability is symbolized by \(\mu_{ij}\), i in the subscript denotes the region, and j denotes the material to be analyzed.

Table 1 Parameter specifications of the material structure
Full size table

3.2 Evaluation, Numerical Solution, and Visualizing of the Magnetic Structure

From the Maxwell equations, a magnetostatic problem is expressed as both vector and scalar magnetic potential “\(A\)” in 2D (\(\left\{ {x,y} \right\} \in Domain\)), noting that \(\nabla .B = 0\). Because the vector potential \(A\) has only one component in 2D, which is in the \(z\) direction \(A_z\) perpendicular to the cross section being studied, the equation is obtained [12, 13] as follows:

$$\left. \begin{gathered} \nabla \times \,(\nabla \times \,{\bf{A}}_z ) = \,\mu_0 \mu_r \,{\bf{J}}_z \hfill \\ \nabla^2 \,(A_z ) = \,\mu_0 \mu_r \,J_z \hfill \\ \end{gathered} \right\}$$
(1)

The current density value \((J_z )\) is calculated from coil current (I), and the cross section of the coil (S), according to the relation between the current, turns of the coil, and the cross section of the coil \((J_z = N\,I/S)\). A uniform current density is perpendicular to the 2D cross section and the vector direction is the same with \(A_z\). Because coding is implemented here, one can adjust the desired number of turns, current, geometry, and ferromagnetic material properties of the system to see how they affect the magnetic vector potential, the magnetic field density, and the magnetic field.

We frequently base our analysis on the magnetic scalar potential, which is related to the mathematical properties of the Laplace equation in 2D (x–y plane). The following scalar equation expresses the solution to the electromagnetic potential problem in the entire region. It is a very convenient equation, particularly studying both the temperature and magnetic effects of the system [14,15,−16] as follows:

$$\nabla .(\nabla A_z )\, = \mu_0 \mu_r \,J_z$$
(2)

Two Dirichlet boundary conditions can be of the following types [11]. We impose conditions (3) and (4) on the air boundary around the solid (just fix potential at one location), as shown in the program scripts, with bc1 and bc2.

$${\bf{B}}_\parallel = 0$$
(3)
$${\bf{A}} = 0$$
(4)

Applying boundary conditions to Maxwell Eqs. (1) and (2) yields identical solutions, as depicted in Fig. 2. The results shown in Fig. 2a and b are produced by executing the "sol1" command, while the outcomes in Fig. 2c and d are generated using the "sol2" command. In Fig. 2a and c, the solution pertains to the solid (magnetic material), while Fig. 2b and d represents the entire domain (whole region). The formulation of a 2D static Maxwell magnetic field is defined through the coding of Pde1, Pde2, sol1, and sol2. Pde1 and Pde2 establish the fundamental structure of the numerical equations, while sol1 and sol2 encapsulate the principal equations along with the corresponding boundary conditions in the x–y planes. The codes required to construct the physical region shown in Fig. 1, a challenging and central aspect of the coding, are included in Appendix. The coding intentionally encompasses two types of materials and their domain behaviors, aiming to encourage students to tackle more intricate problems. The provided codes are ready for execution, allowing one to either utilize them as is or modify them to explore various scenarios. Moreover, even for dynamic nonlinear problems, students can evaluate the impacts of different structural configurations by modifying the codes within the domain.

Fig. 2
figure 2

(a) and (b) Sol1 for the solid and domain regions, (c) and (d) Sol2 for the solid and domain regions for vector potential (\(A_z\))

Full size image
figure a

In the scripts provided above, various combinations and mathematical expressions have been utilized, such as the use of "curl" and derivatives, to define the magnetic field B (both B and B1). This enables a comparison of the results obtained. The term “mu \((\mu = \mu_0 \,\mu_r )\)” and “mu0 \((\mu_0 )\)” are defined in Appendix.

The code snippets are presented here in their original Mathematica forms, aiming to facilitate the implementation of formulas for users. Furthermore, as shown in Table 2, the computational time for both FEM and FDM solutions, under different boundary conditions, indicates only marginal differences. Notably, the Pde1 solution with bc2 yields the shortest computation time.

Table 2 Time used by the different combinations
Full size table

In this study, as illustrated in the code provided above, both FEM and FDM have been employed and are visualized in Fig. 2 (labeled as "sol1" and "sol2," respectively). If the method is not explicitly specified, Mathematica will default to solving the problem using the FDM algorithm. The plotting tools available in Mathematica are just as powerful as the solution methods themselves. Some mathematical requests can be seamlessly integrated into plotting activities. For instance, to calculate vector values on a grid, one can create the grid using the "CoordinateBoundingBoxArray" function, utilizing the following syntax [17]:

figure b

In the analysis above, the values of \(\mu_{r11}\) and \(\mu_{r21}\) have been used for presentation in detail (see Table 1), and the other permeability values are presented in Appendix Fig. 6a, b, and c. In Appendix, Fig. 6a, vector potential is visualized in 3D to show that knowing plotting specialties creates more detail in a numerical solution presentation. The solutions also include some theoretical knowledge. The solutions are identical in either case (sol1 or sol2) inside the solid and the domain. Moreover, the application of boundary conditions (3) and (4) to the solution region satisfies the same solution, and the functions (e.g., vector potential, magnetic field density, and magnetic field intensity), constraints (e.g., boundary conditions), and geometrical entities (e.g., regions, surface) are associated with the mesh of the region geometry.

3.3 Discussion and Results from Comparing Finite Elements Magnetic Method Simulations Solutions

In this section, the results of the values are compared with the FEMM simulation on the ferromagnetic structure. The same dimensions in centimeters shown in Fig. 1 and the same materials shown in Table 1 have been entered into the FEMM magnetic problem analysis, simulation program. The FEMM simulation analysis results in Fig. 3a and b give an analogy with the Mathematica coding program shown in Fig. 4a and b. The relative error of B and H has been defined as approximately 3.1 and 7.2%, respectively, at the same \(\,\,x\,and\,y\) points. Although the values appear to be the same as values, the visualization result, which particularly changes in \(H\) values on the ferromagnetic material structure in Mathematica outputs, is more precise (Fig. 4b).

Fig. 3
figure 3

(a) B values and (b) H values of the solid using FEMM

Full size image
Fig. 4
figure 4

(a) \(B\) values and (b)\(H\) values of the solid using Mathematica

Full size image

In this study, FEM is used in script sol1, and FDM is used in script sol2 to solve the numerical equations. However, using Mathematica tools, the solution method could be easily changed and vice versa or to any other available equivalent numerical or even analytical methods. In addition, visualization of equations gives some advantages to grasping the subject easier. The Mathematica problem statement is given as follows:

  1. 1.

    Initialization: constructing the codes for the physical structure (Fig. 1) to set in the \(x\,,y\) plane (see Appendix).

    • \((\left\{ {x,y} \right\} \in \,\,Domain\,,\,\left\{ {x,y} \right\} \in \,\,solid\,,\,\left\{ {x,y} \right\} \in Coil)\)

  2. 2.

    Solution of Maxwell Eqs. (1) and (2) for the physical structure with the help of Mathematica solver tools (Pde1 and Pd2) to evaluate \(A_z\)(vector potential). In this step, different Maxwell equations may be matched with different boundaries to see the effect of the equations and boundary conditions. Sol1 and sol2 scripts are the solution of the physical problem with boundary conditions (3) and (4). Additionally, H(magnetic field intensity) and B(magnetic field density) can be performed easily in this step, by using the script of B and H in the text.

  3. 3.

    To visualize the problem statement Mathematica has rich tools. A special script is shown in the text. The plotting activities and specialties are shown in the text emphasizing “CoordinateBoundingBoxArray.” The ordinary tool for visualization is shown in Appendix named “StreamDensityPlot.”

The whole numerical analysis processes of the physical problem are illustrated in Figs. 2 and 4. The B and H directions also indicate the current flow direction in the winding (right-hand rule).

Figures 3a, b and 4a, b show the results for \(\mu_{r11} \,,\mu_{r21} \,\) defined in Table 1. The other ferromagnetic materials, such as \(\mu_{r12} \,,\mu_{r22} \,\) and \(\mu_{r13} \,,\mu_{r23} \,\) , are considered with maximum values in Fig. 5a and b with Mathematica outputs and also \(\mu_{r13} ,\mu_{r23}\), as shown in Appendix, Fig. 6a, b, and c. Although the results show high verification, as shown in Fig. 4b, the differentiating ferromagnetic material is more precise in Mathematica for \(H\) values within the region circled as \(\mu_{r21} \,\).

Fig. 5
figure 5

Comparison of (a) B values and (b) H values for different magnetic materials

Full size image
Fig. 6
figure 6

In (a, b), and (c),\(A_z ,B\,,\,and\,\,H\) are figured out for second permeability values defined in Table 1

Full size image

Future research must consider extending models to offer additional features in numerical equations regarding complex geometry structures, such as the following: (1) coupling of potentials and fields, (2) nonlinear formulations involving physical functions with any number of arguments, and (3) application of mixed boundary elements and integral formulations to sophisticated material structures. Additionally, various extensions of the computational tools can be coupled together, which leads to numerous couplings of physical problems and numerical methods (e.g., electric circuit coupling, magnetic circuit coupling, finite element, and boundary element coupling, and time analyses coupling) [12].

4 Conclusion

This educational study provides a detailed and unambiguous solution to the 2D static Maxwell elliptic numerical equation for a specific ferromagnetic material, utilizing a developed program. The study explains the step-by-step solution methods for the Maxwell equations, emphasizing the importance of coding in electromagnetic theory rather than relying solely on simulation programs. The explanation emphasizes the significance of coding, boundary conditions, and numerical solutions in the context of static magnetic Maxwell equations.

Furthermore, the study highlights the following analysis results:

  • The developed program effectively demonstrates the use of two different forms of numerical equations.

  • The solutions within the domain and the solid regions yield similar results, indicating that the outer region of the solid has minimal impact on the solution.

  • Different Dirichlet boundary conditions, represented as A and B, have similar effects on the solution region.

  • The solution time remains consistent across different combinations of equation solutions.

The study is presented in a step-by-step manner to motivate students and researchers to tackle more complex structures and problems. The provided coding is open for further improvement.