## The Boundary Value Problem for Elliptic Equation in the Corner Domain in the Numerical Simulation of Magnetic Systems

**Authors:**Perepelkin E.E., Polyakova R.V., Kovalenko A.D., Sysoev P.N., Sadovnikova M.B., Tarelkin A.A., Yudin I.P.**Issue:**Vol 25, No 3 (2017)**Pages:**253-265**Section:**Mathematical Modeling**URL:**http://journals.rudn.ru/miph/article/view/16210**DOI:**http://dx.doi.org/10.22363/2312-9735-2017-25-3-253-265

#### Abstract

Modern accelerator systems and detectors contain magnetic systems of complex geometrical configuration. Design and optimization of the magnetic systems demands solving a nonlinear boundary-value problem of magnetostatic. The region in which the boundary-value problem is solved, consists of two sub-domains: a domain of vacuum and a domain of ferromagnetic. In view of the complex geometrical configuration of magnetic systems, the ferromagnetic/vacuum boundary can be nonsmooth, i.e. it contains a corner point near of which the boundary is formed by two smooth curves crossed in a corner point at some angle. Thereby, the solution of such a problem has to be found by numerical methods, a question arises about the behavior of the boundary value problem solution around the angular point of the ferromagnetic. This work shows that if the magnetic permeability function meets certain requirements, the corresponding solution of the boundary value problem will have a limited gradient. In this paper, an upper estimate of maximum possible growth of the magnetic field in the corner domain is given. In terms of this estimate, a method of condensing the differential mesh near the corner domain is proposed. This work represents an algorithm of constructing an adaptive mesh in the domain with a boundary corner point of ferromagnetic taking into account the character of behavior of the solution of the boundary value problem. An example of calculating a model problem in the domain containing a corner point is given.

1. Introduction Many physics research facilities use magnetic systems of various configurations, e.g., a system of spectrometric magnets. It is very important to know with a good accuracy the distribution of the magnetic field generated by this system. The problem is actually reduced to formulation of a magnetostatics problem of finding the distribution of the magnetic field generated by the magnetic system under consideration. Since the magnetic system has a complicated configuration, the solution of the problem is usually sought using numerical methods. The domain in which the boundary value problem is solved during calculations of a particular magnetic system often has a piecewise smooth boundary. In this case, the solution of the problem or the derivative solutions can have a singularity. Therefore, the numerical search for the solution requires the use of special methods. 2. Formulation of the Boundary Value Problem The problem to be formulated is the magnetostatics problem of the magnetic field distribution in the corner domain of a ferromagnetic (see Fig. 1). From the Maxwell’s equations and boundary relations (no currents are supposed to be in the region under consideration) it follows that where Ω is the ferromagnetic and the vacuum region (on Figs. 1, 2 region Ω1 = Ω and region Ω2=Ω ); Γ is the boundary; and ⃗ strength vectors, respectively. and ⃗ are the magnetic field induction and Figure 1. The corner domain Figure 2. The angular sector For the ferromagnetic region Ω2, we can write ⃗ = 0()⃗ , where = ⃗ , () is the permeability, and 0 is the vacuum permeability. For the vacuum region Ω1, we have ⃗ = 0⃗ . Since there are no current sources in the region Ω = Ω1 Ω2, the field is potential, and thus the following representation holds: where () is the scalar potential. The consequent formulation of the boundary value problem is where the function () satisfies the conditions: 1. ; 2.; 3. Let us consider the function ¯(), an analogue of the function (), for which conditions 2 and 3 are replaced by ¯ (′) = 1 for ′ ): 0, where 0 is “large enough”. In what follows we will assume that the solution to (1) is ∈ (Ω ∪ Γ ∪ Γ1 ∪ Γ2), and thus it follows that ∃0 > 0 ∀ ∈ Ω ∪ Γ ∪ Γ1 ∪ Γ2 : |()| < 0. 3. On a Certain Boundary Value Problem Before proceeding to the main statements of the paper, we consider an auxiliary problem that is discussed in detail in [1], namely, the boundary value problem (see where Ψ ∈ (1) (Γ), = 1, 2, (see Fig. 2). Let us introduce the polar coordinate system. The solution is sought by the method of separation of variables: ∼ () · Φ(). As a result, we have Φ′′ + 2Φ = 0, 2′′ + ′ - 2 = 0. Thus, by virtue of boundedness of () at the origin, the solution for () will be () , and for Φ() there will be eigenfunctions divided into two groups, symmetric about = 0 and antisymmetric about = 0. In the former case, the eigenfunctions take the form Here the constant is determined from the boundary ratio for normal derivatives and the eigenvalues can be determined using the continuity relation for the solution (, ) at the boundary Thus, either tg = 0 and hence = 4, or = smallest root of the equation ± 4, where 1 is the (2) The singularity is introduced in the solution by the series term 1 Φ From (2) it follows that . (3) This means that if 1 = 2, |()| will be limited. 4. Behavior of the Solution of the Boundary Value Problem Let us consider boundary value problem (1) with the permeability function ¯, region Ω (see Fig. 1). Statement 1. ∃ > 0 ∀ > 0 0 < (, ) < : |()| < , where (, ) is the distance between the points and , and by boundedness on Γ is meant boundedness on Γ+ and Γ-. Proof. We will prove it by contradiction. Let us assume that it is not true. Then ∀ > 0 ∃ > 0 0 < (, ) < : |()| ): . Let us take = max (0, 40√), 0 < 0 < , then for : 0 < (, ) � 0 the condition should hold. |()| ): 0 ⇒ ¯ () = 1. (4) We introduce a polar coordinate system with the origin at the point . Let (, ) be the solution of our boundary value problem satisfying condition (4). Then on Γ it should satisfy the conditions , (5) and by virtue of continuity of (0, ) (0, )|=0+ = (0, )|=0- = (0, 0) (6) should hold. It follows from (5) and (6) that Thus, putting in the 0-vicinity of the point , we obtain the boundary value problem , (7) where Ψ ∈ (1) (Γ). From (3) (and also from [1]) we find that (7) has no singularities, i.e., lim |()| = √︀2 + 2 � 40√ = , where 1 and 1 are the Fourier series coefficients for the function () at the boundary of Γ. Consequently, we have arrived at a contradiction, which proves our statement. Thus, it follows from the statement that the magnetic field is bounded in the corner domain provided that the permeability function satisfies the conditions: 1) ¯ () ∈ (1)[0, +∞); 2) ∃0 > 0 ∀′ ): 0 : ¯ (′) = 1. Note an interesting fact. Let us solve boundary value problem (1) and let its solution have unbounded . This means that in the vicinity of the point the permeability function ( ) will tend to unity. Since the number of figures in the mantissa is limited, it will turn out that in a certain small vicinity of the point the function ( ) will be equal to 1. That is, there arises boundary value problem (1) with the permeability function ¯ () that has bounded and thus we get a contradiction. Consequently, numerical calculations cannot “theoretically” yield a solution with the infinitely growing , and we will seek the solution of another boundary value problem, namely, problem (1) with ¯ (), where is limited. But the solution of problem (1) with ¯ () does not coincide in the general case with the solution of initial problem (7) with (). It is therefore necessary to use special methods for solving this problem. One of these methods is considered in [2, 3] for the solution of the equation div [ ( ) ] = 0 in the corner domain. 5. Estimation of the Magnetic Field Growth Let us show that the magnetic field in the corner domain of a ferromagnetic satisfies the condition , (8) where 0 is a constant; () is a bounded function; and is the distance to the corner. The integral formulation of the magnetostatic problem allows the magnetic field to be represented as , (9) where ⃗ is the field from the current sources; ⃗ is the ferromagnetic magnetization 1 vector; the function Ψ (, ) is equal to or ln sp for the three-dimensional 2 and the two-dimensional case, respectively; and Ω is the ferromagnetic domain (see Fig. 2: region Ω1 = Ω and region Ω2 = Ω ). The magnetization vector is defined as ⃗ = 0 () ⃗ = 0 ( () - 1) ⃗ , where 0 is a constant; () is the magnetic susceptibility; and () is the permeability of the ferromagnetic. Given high fields ( → ∞), the representation () = 1 + - 2 , → ∞ is valid, where and are positive constants. Consequently, when → ∞, = ⃗ is limited by a constant 0 = 0. Let us consider the 2D case. From (9) we obtain Here the first term is limited, and we therefore estimate the second term We calculate the integral where = Ω ∩ () is the angular sector at the corner point (see Fig. 2). The integral over the domain Ω / will be limited, and we therefore consider only the integral over the domain Then we use the expression for the generating function and obtain where 1 is a constant, and 1 () is a bounded function. Thus, the validity of expression (8) is ascertained. 6. Method of Mesh Condensing in the Corner Domain In [4,5] there are examples of constructing a differential mesh for some boundary value problems in corner domains. The main idea is to condense the differential mesh or finite elements for obtaining an admissible problem approximation error. This error involves integrals over elementary domains estimated by the quantities of the form ‖‖ , where is the diameter of the elementary domain or mesh cell; is a positive number; ‖‖, is the norm of the function with the th derivative in this domain; and is a constant independent of all these factors. Then we can require, for example, that quantities ‖‖ be identical in the domain under consideration. To this end, can be decreased in inverse proportion to ‖‖, on approach to the singular points. We demonstrate the validity of the following statement. Let ⃗ () be the solution of the magnetostatic problem in the integral formulation found by a numerical method and estimate is valid: ⃗ () be the exact solution. Then the following where 1, 2, and 3 are positive constants and is the diameter of the domain , which is a differential grid cell containing the ferromagnetic corner. By virtue of (9), the following expression for ⃗ () holds: where is the field in the cell Ω , = 1 . . . ; Ω = Ω ; and is the distance from the point to the point ∈ Ω . We consider the difference Since the quantity ⃗ < 0 is limited, it follows that ⃗ ⃗ - ⃗ ⃗ ( ) < As a result, using the estimate obtained above, we arrive at the expression It remains to estimate ⃗ - ⃗ , where the domain is the of the corner point . Using (12), we obtain () - -domain of the corner point . Using (12), we obtain where = 2, and 1, 2, and 3 are constants. The result corresponds to (10). Based on the aforesaid and inequality (10), we propose a differential mesh condensing method Here 0 is a constant; is the number of partitions along the coordinate axis ( or ) in the corner domain; is the grid spacing; and is the coordinate of the grid node along the or axis (the origin of the coordinates is at the corner point), || < 1. 7. Calculation of the Magnetic Field in the Corner Domain In Section 5 we gave the upper estimate of the admissible magnetic field growth in the corner domain of a ferromagnetic. In Section 6, based on this estimate, we proposed a method of condensing the differential mesh in the corner domain. In this section we present numerical calculations of a magnetic system using this method. It is evident from the calculations that the above mesh condensing method substantially improves the accuracy of the calculated magnetic field distribution. We considered a magnetic system depicted in Fig. 2. The domain Ω corresponds to the vacuum, and the domain Ω is filled with a ferromagnetic. A boundary value problem corresponding to the formulation of the magnetostatic problem with respect to the vector potential , was solved in the domain under consideration. The function is defined as = 1/ (), where () is the permeability of the ferromagnetic; is the modulus of the magnetic induction vector ⃗ = rot ⃗; Γ± is the interface; and is a constant. The efficiency of the differential mesh condensing method described in Section 6 was estimated by the following calculations: Variant 1. The solution to problem (13) was found on a sequence of (10) (10) meshes, where = 1, 2, 4, 8, 10, 20, 40; that is, 10 10, 20 20, . . . , 400 400 meshes were obtained. The mesh spacing in the corner domain Π was uniform. Variant 2. Problem (13) was calculated on the same sequence of meshes at = 1, 2, 4, 8, but the mesh spacing in the domain Π was chosen using the differential mesh condensing method described in Section 6. In the domain Ω/Π the mesh spacing was not changed as compared with the previous variant. The number of node points in the domain Π was the same, and only their distribution was changed. The results of the variant 1 calculations were taken to be a sort of reference because the accuracy of the calculated solution was assumed to increase with increasing number of partitions, except probably for a particular corner domain. Then the results of the variant 2 calculations were compared with the results of the variant 1 calculations. Figures 3-6 shows distributions () at = 1 for different meshes. In all figures the reference distribution () calculated on the 400 × 400 mesh is designated as trace 1. The plots of traces 2 and 3 are the distributions () calculated by variants 1 and 2, respectively. The distributions in Figs. 3-6 are calculated on the 10 10, 20 20, 400 400, and 80 80 meshes, respectively. It is evident from Figs. 3 that the accuracy of the variant 2 calculations (nonuniform mesh) is substantially higher than that of the variant 1 calculations (uniform mesh). Thus, it follows that the proposed method of constructing a differential mesh in the corner domain is worth using and yields results comparable in accuracy with the results obtained only on meshes with the number of nodes in each axis four to five times greater than in the initial mesh.An algorithm of thickening differential mesh near the corner point has been developed. It allows one to significantly reduce the computation time and simultaneously to increase the accuracy of the solution of the boundary value problem. Variant 1 - maximum of relative error is 11.085%; variant 2 - maximum of relative error is 1.091%. Figure 3. Distributions () at = 1. Mesh 10x10 Figure 4. Distributions () at = 1. Mesh 20x20 Figure 5. Distributions () at = 1. Mesh 400x400 Figure 6. Distributions () at = 1. Mesh 80x80 8. Results of Modeling of Some Magnetic Systems The significance of numerical modeling at the investigating of magnetic systems is defined by not only known dignities of computational experiment but also that the measurement of magnetic field is labor-intensive and expensive problem. The process of the mathematical modeling of magnetic systems (see Figs. 7 and 8), as the authors of this work have presented, should be divided into two large stages. Figure 7. General view of spectrometric magnet 1SP-40-4B Figure 8. Type of solenoidal magnet 1. Results of modeling magnetic systems SP-94 and 1SP-40-4B In the experiment performed at the Veksler and Baldin Laboratory of High Energy Physics, JINR, the SP-94 [6] magnet is used. It was necessary to select the configuration of the core and the current coils such that the quantity (0, 0, ) d had the maximum value. Fig. 9a presents the distribution of in -∞ plane for this configuration of the magnet = 94 (variant 1). First, the initial configuration (variant 1) was calculated. Here (0, 0, ) d = 2.314, where = 1.5 m is the dimension along the axis of the region where the magnetic field was calculated for variant 1. Second, the initial configuration (variant 2) was calculated. Fig. 9b presents the distribution of in the plane for configuration of the magnet = 94 (variant 2). Here (0, 0, ) d = 2.987, which is 1.291 times greater than for the initial configuration in variant 1. Figure 9. Distribution of for three configurations of the magnet SP-94 Two configurations of the magnet 1SP-40-4B NIS (variant 1, 2, respectively) for which the numerical computations of the magnetic fields were performed, i.e., in fact, a nonlinear inverse problem of magnetostatics were solved. The purpose of the simulation is to find by a calculation method such a magnet geometry that the region of the homogeneity of the magnetic field would be essentially larger as compared to the existing magnet configuration. Fig. 10 present distributions of the components of the magnetic field for our two variants of configurations of the magnet 1SP-40-4B in a 3D case. Clearly, for variant 2 the distributions of the components got more smooth. Figure 10. Dependences (, 0.3, ) for two configurations of the magnet, current 1100 A 2. The Solenoid Type Magnetic Field Detector Modeling Magnetic systems are very important parts [7, 8]. To create the necessary configuration of magnetic field, the repeated solution of nonlinear boundary value problem of magnetostatics is needed. In the present work, we consider the problem of creation of homogeneous map of magnetic system of solenoidal type (see Fig. 8). As a result of optimization, the geometric parameters of magnetic system were chosen in such a way so as to get maximal size of the domain of homogeneity of the magnetic field. Due to symmetry, only 1/24 part of the geometry with corresponding boundary conditions is modeled. The calculations were performed (using two software products: TOSCA and native MFC) by the method of finite elements on tetrahedral mesh with 5 000 000 elements. In Figs. 11 and 12 the domains with the degrees of homogeneity of magnetic field are of 0.1% and 0.5%, respectively. The black continuous line shows that the homogeneity of 0.1% is needed. In Fig. 11, the scale of magnetic field has site from 0.99-1.001 T, in Fig. 12 from 0.998-1.002 T. Figure 11. Field homogeneity is ±0.1% Figure 12. Field homogeneity is ±0.5% 9. Conclusions The upper estimate for the admissible growth of the magnetic field ⃗ () in the corner domain where 0 is a constant, () is a bounded function, and is the distance to the corner, is asymptotically obtained for the case of () → 1 when → ∞. A method of condensing the differential mesh in the corner domain is proposed, which appreciably improves the accuracy of the calculated solution. The numerical modeling results are presented for the SP-94 magnet system in the Delta-Sigma experiment performed within the Topical Plan of JINR on international collaboration. Two-dimensional and three-dimensional modeling of the configuration of the magnet core and current coils was performed to obtain the maximum value of the integral By a numerical method a configuration of the magnet 1SP-40-4B VBLHEP, JINR has been selected for which the width of the domain of the homogeneity of the magnetic field has grown up from 0.5 m to 1.0 m, i.e., twice. This growth considerably increases the accuracy of regenerating the pulses of decay particles in the physical reaction under study (search for pentaquarks). As a result of optimization, the geometric parameters of the solenoid type magnetic field detector were chosen in such a way so as to get maximal size of the domain of homogeneity of the magnetic field.

### E E Perepelkin

Lomonosov Moscow State University
**Author for correspondence.**

Email: perepelkin.evgeny@phys.msu.ru

GSP-1, Leninskie Gory, Moscow, 119991, Russian Federation

### R V Polyakova

Joint Institute for Nuclear Research
Email: polykovarv@mail.ru

6 Joliot-Curie St., Dubna, Moscow region, 141980, Russian Federation

### A D Kovalenko

Joint Institute for Nuclear Research
Email: kovalen@dubna.ru

6 Joliot-Curie St., Dubna, Moscow region, 141980, Russian Federation

### P N Sysoev

Lomonosov Moscow State University
Email: apc_box@mail.ru

GSP-1, Leninskie Gory, Moscow, 119991, Russian Federation

### M B Sadovnikova

Lomonosov Moscow State University
Email: apc_box@mail.ru

GSP-1, Leninskie Gory, Moscow, 119991, Russian Federation

### A A Tarelkin

Lomonosov Moscow State University
Email: tarelkin.aleksandr@physics.msu.ru

GSP-1, Leninskie Gory, Moscow, 119991, Russian Federation

### I P Yudin

Joint Institute for Nuclear Research
Email: yudin@jinr.ru

6 Joliot-Curie St., Dubna, Moscow region, 141980, Russian Federation

- G. Strang, G. Fix, An Analysis of the Finite Element Method. Second edition, Wellesley-Cambridge Press, 2008.
- E.P. Zhidkov, E.E. Perepelkin, An analytical approach for Quasi-Linear Equation in Secondary Order, Computational Methods in Applied Mathematics 1 (2001) 285–297.
- E.E. Perepelkin, R.V. Polyakova, I. P. Yudin, The Boundary Value Problem for Elliptic Equation in the Corner Domain, Bulletin of Peoples’ Friendship University of Russia (2) (2014) 410–414, in Russian.
- E.A. Volkov, Method of meshes and infinite domains with a piecewise smooth boundary, Dokl. Akad. Nauk SSSR 168(3) (1966) 978–981, in Russian.
- V.V. Shaidurov, Numerical Solution of the Dirichlet Problem in a Domain with Angles, Nauka, Novosobirsk, 1982, in Russian.
- P. Yudin, V.A. Panacik, E.E. Perepelkin, R.V. Polyakova, A.N. Petersky, Peculiar Features of Numerical Modeling of the Modified Spectrometer Magnet Field, Computer Research and Modeling 7 (1) (2015) 93–105.
- E. Perepelkin, et al., The ATLAS Experiment at the CERN Large Hadron Collider, Vol. 3, Aad, JINST, 2008.
- E. Perepelkin, et al., Commissioning of the Magnetic Field in the ATLAS Spectrometer, Vol. 177–178, 2008.