CC BY
28BOUNDARY ELEMENT METHOD IN 3D PROBLEMS OF MATHEMATICAL PHYSICS
V. Ivanov
Institution of Computational Technologies SB RAS, 630090, Novosibirsk
UDC 519.6
DOI: 10.24411/9999-016A-2019-10027
Numerical algorithms of the boundary problems described by Laplace, Poisson, Helmholtz and Maxwell equations based on the integral representations are described. Combinations of single-layer, double-layer and volume potentials are used to represent the solution. Analytical technique of singularity extraction is used to obtain precision and stable numerical solutions for a potential and its high-order derivatives. The numerous examples of test problems and applications to the precision problems of electron optics, accelerator physics and plasma-beam interaction are demonstrated. The original algorithms can be used for the problems of analysis, optimization and synthesis of physical electronic devices. Decomposition algorithms to solve complex three-dimensional problems are presented. keywords: Boundary Element Method, integral representation, single-layer potential, double-layer potential, numerical approximation, singularity.
Introduction
The mathematical apparatus for describing natural phenomena is usually called the equations of mathematical physics. These equations admit different forms of representation. The most well-known are writing in the form of partial differential equations, variation formulations of problems and integral representations. Thus, the stationary problems of electromagnetism, hydrodynamics, and thermal physics lead to boundary value problems for the Poisson equation
= p, (1)
where A — is the Laplace operator, $ — is a potential, but p — is a density of field sources. The time-dependent equations of electrodynamics are described by the system of Maxwell's equations
dB ~dt dD
— -r otH = J, (2)
dt w
Here E and H — are the strength of electric and magnetic fields, D = eE and B = ^H — are the vectors of electric and magnetic induction, J — is a current density of excitation field, c = — is a speed of light in vacuum, e and ^ — are the permittivity and permeability of a physical environment. To obtain the integral representation, we use the second Green's identity
\dS.
+ r ot E = 0,
— r ot H = J,
divB = 0,
di vD = 0,
JJJ (uAv - vAu)dV = jj) [u^- - v^ dS, (3)
which is valid for any two harmonic functions u and v. Here, on the right-hand side of the identity, integration is over the boundary S of the volume V occupied by the field, and n is the internal normal vector to this boundary.
ISBN 978-5-901548-42-4
If as the trial function u we choose the potential $(p) at some arbitrary point p, and put the function v equal to v = 1/Rpq, where Rpq — is the distance between the observation point p and the source point q, over which the integration in formula (3), then, taking into account the relation A-J— = S and the property of the ¿-function, we obtain the required representation
*(p) = $ M^ihdSq + ^dVq. (4)
_±_dS -£ + [ff^c
' dnq Rpq q JJ dnq Rpq q JJJ Rpq
In this representation, the first integral is called the double-layer potential, the second is the simple-layer potential, and the third is called the Newtonian potential.
The integral representations for the system of Maxwell equations have the form
j>M = -/
BdS = 0, (5)
HH = JJ + J\ dS, DdS = iff pdV.
The properties of integral equations was represented in fundamental monographs by C. Baker [1]. In the numerical solution of the equation of mathematical physics, they are approximated by some set of orthogonal functions Hm. So that the norm of the difference of the original functions H and the approximated ones satisfies the relation ||H — Hm\\ « o(hm), where h — is the characteristic scale of the discretization of the equations, and the number m is called the order of approximation. After sampling, the original equations reduce to a system of linear algebraic equations.
The use of non structured mesh with some geometrical elements on the boundary surface produces the Boundary Element Method (BEM) The technology of development of numerical algorithms of BEM is considered in monographs [2]- [4], and the convergence of these algorithms is investigated in [5]. Consider a mixed boundary-value problem for the Poisson equation (1). Let the Dirichlet boundary condition be given on the part Si of the boundary of the domain
$ |si = Uo(p) ,p eSi, (6)
but at the rest part of the boundary S2 the Neumann condition is given
— g S =Eo(p) ,p eS2. (7)
Here U0 and E0 — are given continuous functions, and n is the inner normal to the boundary. Applying the boundary condition (6) to the integral representation (4), one can obtain the Fredholm integral equation of the first kind with respect to the unknown function a, but of the second kind equation with respect to the unknown function ^ at the boundary
ii ^^ l^dSq + 2^(p)+/7 a^dSq = fi(p),p eSi, (8)
JJs2 anq Rpq JJs 1 Rpq
if d d 1 ff d 1
£m, an: TU,dS' + 2™<p> +1 a<" an, RZ,dS< = «p>'p e <9>
where
h(p) = U (p) - JJJ^^dVq ,P eSi, (10)
h(p) = -Eo(p) -JJjvP(i)JTdVq,P e S2. (11)
The kernels of the integral equationd have the singularities at p ^ q . The regularization procedure consists in adding a term to the integrand, the integral of which can be calculated analytically. In this case, the equation in the neighborhood of the singularity So, where p ^ q becomes
ff dsq + a(p} ff ldsq = w (12)
JJS 0 Rpq J J S0 Rpq
The integrand in the first integral has no singularity, and the second integral can be calculated analytically when the fragment of the surface S0 is represented by a second-order surface. Thus, equation (12) is an integral equation of the third kind.
1 Three-dimensional algorithms
When the BEM implementing for three-dimensional problems, a dilemma arises, using structured or non-structured meshes when discretizing the boundaries of the region. We demonstrate the effectiveness of numerical algorithms for both versions [6]- [7]. Let the boundary of the domain consist of m smooth surfaces S = |J™ 1 S\ , each of which admits a parametric representation
v = Vi(T, r/),y = yi(t, r/),z = Zi(t, r/),T € [ai,fti],r/ € [71, ¿¿]. (13)
Details of algorithms for solving three-dimensional problems are described in the monographs [6-7]. So in the case of a bilinear spline for surface sources
a(T,V) = {[ai+i,j+i(T - Ti)+ai,j+1(Ti+1 - t)](i1 - r]j) + [ai+hj(t - ri) + aiiJ-(Ti+1 - T)](ilj+1 - i1)} /(hThv),
(14)
where hT = Ti+1 - Ti, hv = rjj+1 - r/j finally we get the linear algebraic system GX = F with the global matrix G, unknown vector X and right-hand vector F, which includes the boundary conditions. Matrix elements are given with the formula
G« = i^tj ' <15>
Here — are the moments of the bilinear form in the approximation of the field sources, J — Jacobian of the coordimate transformation.
Let us consider the non uniform parallelepiped mesh in Cartesian coordinates {xi} x {yj} x {Zk} with linear approximation for the space charge [9]- [10]. In order to evaluate the potential:
^(x,У, Z) = ^f [ [ -^ -7TdxdVdZ. (16)
Jy, J*k R(x -x ,y - y',z - z )
One need the four moments of space charge
J1 = IK Rdxdydz,J2 = JJJ xR^dxdydz, J3 = JJJ -x^dxdydz, J4 = JJJ Rdxdy,dz, (17)
where r2 = (x - x')2 + (y - y')2 + (z - z')2. The other integrals are obtained by cyclic change of variables x, y, z. All integrals can be evaluated analytically.
The CPU-time T for the boundary problem solving can be represented with the formula T = aN3 + bN2, where N — is the dimension of matrix G, but a and b — are some numerical constants. First term of this formula is a time for solving of the linear system by the Gauss elimination method, and second term is a time for evaluation of the matrix elements. First term dominates for large enough N. To overcome this strong cubic dependence, H.Schwarts suggested the decomposition idea, when a large problem is spitted onto K, overlapped blocks separated each other with the auxiliary boundaries. Then K boundary problems of lower dimensions are solving in the iterative loop, where the boundary condition on the auxiliary bounds are evaluated from the solution of neighbor block at previous iteration
We can make an estimation for the efficiency of decomposition algorithm comparing with the direct solution of the whole problem of dimension N. Let in simplest case all blocks have equal dimensions, i.e. N = Km. Then the relative efficiency for large N can be estimated by the formula
E = aN3 + &N2 ^ J2 IK (am3 + bm2) k ,
where — is a total number of iterations to reach the predefined accuracy . Publication [11] shows that typical value I = 3 ^ 4 for e = 0.1%.
Another approach is based on the use of an unstructured surface mesh. In this case, first the boundary surface is transformed into a set of geometric elements in a simple form, as in the FEM. Polynomial approximations are set on the surface of the elements, and the natrix elements are calculated. In many cases, integrals over the surface of elements can be calculated analytically. For example, in the works of M. Okon and R. Harrington [12]- [14] analytical expressions for the potential of triangular elements with piecewise constant and linear distribution of the charge density on the surface of an element.
We also give an example of the solution of the equations of electrodynamics in the case of harmonic oscillations. We write down the integral equations for the component HT, Hr of the magnetic field [8]:
H(to, vo) = j) {Gt0,vHt(t, n) — Gt0,v(t, V)]dS, (19)
Hr(T0,V0) = Q <js {Gno,riHT(T,11) —Gn0,ri(T, 1])}dS.
Here (t0, r/0) e S — are the values of variables t, r/ at the observation point , Q — is the solid angle at which the surface S is visible at this point. The tensor Green's function Ga^(t0, r/0; t, rf) = [oiVfyP, $ = cod(jR)/R, k = w/c is an elementary current component j = oS(r0)exp[«(kr0 — wt)], r0 = (t0,^0) , oriented in the direction of the vector o and varying in harmonic law with a frequency w, c — is a speed of light, R = ||r0 — r\\ — is a distance between the observation point r0 and the integration point r = (t, r/) . The system of integral equations (19) can be written in the operator form [I — L(w)]Ht = 0, where I — is the unit operator, but Ht — is the component of the magnetic field vector tangential to the surface S. With numerical approximation, this equation reduces to a system of linear equations A(w)Ht = 0 with square matrix A, which approximate the operator I — L(w). A nontrivial solution of this system exists if det(A) = 0.
If bilinear approximation is used for an unknown vector-valued function Ht, the boundary collocations method reduces the original problem to a linear system with the matrix elements
Ga,ß —
4 Í
/ (To,Vo; T, r])JdTdq. (20)
The function Ga,ß has an integrable singularity, for the separation of which the integral (19) is represented in the form
/ ^Ga,ßJdrdß — / (Ga,ß — Ga,ßJdrd-q + / ^mGaßJdTdq, (21)
where Gaß — lim^o Ga,ß The integrand in the first integral on the right-hand side does not contain a singularity, and the second integral for quadratic surfaces is calculated analytically.
2 Solving practical problems
As a practical problem, let us consider the design of a gun of an 80-megawatt klystron with a sheet beam [16]-[17]. Parameters of the gun: current 250A, accelerating voltage 415 kV, working frequency of klystron 91 GHz. In traditional klystrons with a cylindrical beam, the effect of the space charge creates serious obstacles to constructing devices with a power of more than 75 megawatts, so more powerful klystrons can be built on new principles: multi-beam or with a sheet beam. In multi-beam klystrons, when bundling individual beams into a single bundle, it is difficult to obtain a beam with a small transverse phase volume, therefore, the production of laminar beams has prospects in a plan with a sheet beam. The geometry of the gun of such a klystron is shown in Figure 1.
Figure 1: Diode gun on left ( 1 — cathode, 2 — focusing electrode, 3 — anode) and triode gun (on right) of powerful sheet-beam klystron
Figure 2: Klystron gun simulation on structured mesh
Figure 3: Klystron gun simulation on unstructured mesh
A significant beam current and a high degree of relativism lead to the problem of a self-consistent field in which it is required to take into account the Coulomb field of the beam charge and the self consistent magnetic field of the relativistic beam.The results of the calculations are presented in Figures 2 and 3. The first one is based on a structured grid. In the calculation shown in Figure 3, the electrode surfaces are approximated by a set of triangles that are milled in regions with small radii of curvature. For a beam with a current of 250 A, it took about 10 iterations to obtain current convergence with an accuracy of 0.5
As an example of a practical problem of calculating three-dimensional electrodynamic structures [8], a cell of a biperiodic structure with an inductive coupling was depicted in Fig. 4. The results of experimental measurements of the resonance frequency on the type of oscillations k/2 corresponded to 430 MHz, and numerous calculations
for a surface partition by 40 elements gave the value of this the value of 434 MHz. When the number of elements was enlarged to 826, the relative error in calculating the natural frequencies was 5 • 10-4
Figure 4: A cell of bi-periodic RF-structure
References
[1] Baker C. T. The Numerical Treatment of Integral Equations. Oxford: Oxford University Press, 1977.
[2] Banerjee P. K. et al., editors. Developments in Boundary Element Methods, Vols 1,2,3, Applied Scien;ce Publishers, London, 1979.
[3] Brebbia C. A, Telles J. C, Wrobel L. C. Boundary Element Techniques. Berlin: Springer-Verlag; 1984.
[4] Mackerle J, Andersson T. Boundary element software in engineering. Advances in Eng. Software. 1984; 6: 66-102.
[5] Wendland W. L. Boundary element methods and their asymptotic convergence. In P. Filippi, editor. Theoretical Acoustics and Numerical Techniques, CISM Courses. Wien, New York: Springer-Verlag; 1983: 135-216.
[6] Ivanov V. CAD methods for the electronic devices. Part 1: methods for calculating physical fields. Novosibirsk: Inst. of Mathematics; 1986. 194 pp. In Russian.
[7] Ivanov V. Computational methods, optimization and synthesis in electron optics. Hmbg: Palmarium Academic Publishing; 2016. 525 pp.
[8] Ivanov V, Teryaev V, Karliner M, Yakovlev V. Application of the method of boundary integral equations for the calculation of high-frequency resonators. Zhurnal Vychislitelnoi Matematiki I Matematicheskoi Fiziki. 1986; 12 :.1900-1905. In Russian.
[9] Ivanov V. Green's Function Technique in Forming of Intensive Beams. Int. J. of Modern Physics A, 2009; 24(5): 869-878.
10] Ivanov V. Analytical Technique in the Boundary Element Method for 3D Problems of Electron Optics. Microscopy & Microanalysis. 2015; 21(4): 236-241.
11] Ivanov V. Decomposition method for complex 3D problems in electron optics. J. of Physical Science & Applications. 2015; 5(2): 96-100.
12] Okon M, Harrington R. F. The potential due to a uniform source distribution over a triangular domain. Int. J. for Numer. Meth. in Engng. 1982; 18: 1401-1419.
13] Okon M, Harrington, R. F. The potential integral for a linear distribution over a triangular domain. Int. J. Numer. Methods Engng. 1982; 18: 1821-1828.
14] Okon M, Potential integrals associated with quadratic distribution in a triangular domain, Int. J. Numer. Methods Engng. 1982; 21(2): 197-209.
[15] Ivanov V. The method of analysis of three-dimensional nonstationary flows of charged particle. In: Numerical analysis. In: Proc. Inst. of Mathematics SB RAS. Novosibirsk, Nauka. 1989; 15: 172-187. In Russian.
[16] Ivanov V, Krasnykh A. Analytical computations in solution of 3D problems of physical electronics by BEM. In: Proceedings of the ACES-2002, March 18-22, 2002, Monterey, CA.
[17] Ivanov V, Krasnykh A. Scheitrum G, Sprehn D, Ives L, Miram G. 3D modeling activity for novel high power electron guns at SLAC. In: Proceedings of the Particle Accelerator Conference (PAC 03), Portland, Oregon, 12-16 May 2003: 3312-3314.
Ivanov Valentin Yakovlevich — Doctor of Sciences (Physics and Mathematics),
Institute of Computational Technologies SB RAS;
e-mail: [email protected]. Received — 30 April 2019.
