A Nonlocal Operator Method for Partial Differential Equations with Application to Electromagnetic Waveguide Problem

Timon Rabczuk , Huilong Ren and Xiaoying Zhuang

Abstract: A novel nonlocal operator theory based on the variational principle is proposed for the solution of partial differential equations. Common differential operators as well as the variational forms are defined within the context of nonlocal operators. The present nonlocal formulation allows the assembling of the tangent stiffness matrix with ease and simplicity, which is necessary for the eigenvalue analysis such as the waveguide problem.The present formulation is applied to solve the differential electromagnetic vector wave equations based on electric fields. The governing equations are converted into nonlocal integral form. An hourglass energy functional is introduced for the elimination of zeroenergy modes. Finally, the proposed method is validated by testing three classical benchmark problems.

Keywords: Nonlocal operator method, Variational principle, Nonlocal operators,Hourglass mode.

1 Introduction

For the analysis of complex structures in engineering, mesh generation remains a laborious and time-consuming process, which often requires human interaction with meshing program and corrections for local mesh. To alleviate the mesh entanglement, socalled meshless or meshfree methods have been proposed during the 1990s [Viana and Mesquita (1999); Xuan, Zeng, Shanker et al. (2004); Razmjoo, Movahhedi and Hakimi(2011); Nicomedes, Bathe, Moreira et al. (2017)]. Meshless methods have been afterwards developed, enriched and applied to analyze a wide variety of problems in engineering including electromagnetic problems. For example, Viana and Mesquita applied the meshless Moving Least Square Reproducing kernel Method to solve twodimensional static electromagnetic problems. Liu et al. [Liu, Yang, Chen et al. (2004)]proposed an improved formulation of the element-free Galerkin method for electromagnetic field computations and studied the selection of the weight function, the treatment of imposing boundary conditions and interface conditions. The present original nonlocal operator method (NOM) proposed by the authors can be used to solve the partial differential equations by replacing the local operators in PDEs with newly defined nonlocal operators. The nonlocal operator is defined in the integral form based on the nonlocal interaction but converges to the local operator in the continuous limit. The nonlocal operator method is consistent with the variational principle and the weighted residual method, based on which the tangent stiffness matrix can be obtained with ease.The nonlocal operator method can be viewed as a generalization of dual-horizon peridynamics [Ren, Zhuang, Cai et al. (2016); Ren, Zhuang and Rabczuk (2017)] or peridynamics [Silling (2000); Silling, Epton, Weckner et al. (2007)]. In this paper, we develop and apply the nonlocal operator method to electromagnetic problem.

Electromagnetic analysis has been an indispensable part of many engineering and scientific study since Maxwell established a unified electromagnetic field theory-the Maxwell equations-in the 19th century. The Maxwell equations describing electromagnetic waves have numerous applications including radar, remote sensing,bioelectromagnetics, wireless communication and optics, just to name a few. Several computational methods have been developed for the solution of the Maxwell equations including the method of moments [Gibson (2007)], finite element method [Jin (2015)],time domain finite difference method [Taflove and Hagness (2005)], ray theory[Deschamps (1972)], meshless/meshfree methods [Ho, Yang, Machado et al. (2001)],asymptotic-expansion methods [Bouche, Molinet and Mittra (2012)] and eigen expansion method [Chew, Jin, Lu et al. (1997)]. The finite element method and time domain finite difference method can capture complex shapes and inhomogeneous materials while the moment method is well known for its high precision and efficiency for mainly linear problems and simple geometries [Jin (2015)]. The Finite-Difference Time-Domain(FDTD) method [Yee (1966)] is a method for directly solving the Maxwell equation in the time domain. The FDTD method has merits such as the explicit time integration, high efficiency in storage and natural parallelization. The method of moments (MoM) or boundary element method (BEM) is a numerical computational method of solving linear partial differential equations which have been formulated as integral equations (i.e., in boundary integral form) [Harrington (1993)]. This method possesses high accuracy but requires artificial intervention to handle the integral equations. In addition, this method is only applicable for regular shapes and homogenous materials. The multi-layer fast multipole technology based on the moment method is often employed for the purpose of computational efficiency. The accuracy of the finite element and FDTD is lower than the moment method since both finite element and time domain finite difference have numerical dispersion errors, which does not occur for the moment method.

The purpose of this paper is to develop a framework of nonlocal operator method exploiting variational principles and to reformulate the electromagnetic governing equations. Therefore, the local differential equation is converted into nonlocal integral form. The content of the paper is outlined as follows: The governing equations for electromagnetic fields in the time domain and frequency domain are described in Section 2. The concept of nonlocal operator method including definitions of the nonlocal curl and gradient operators are introduced in Section 3. Furthermore, the nonlocal formulation based on the variation of nonlocal operators in discrete form are presented in details to finally obtain the consistent tangent stiffness. In Section 4, we convert the differential equation into the nonlocal form. In Section 5, the hourglass energy functional to remove zero-energy modes is proposed, which are inherited from particle-based formulations based on nodal integration. The residual and tangent stiffness matrix of the hourglass functional is also derived. The nonlocal integral form of the electromagnetic equations in the time-domain is derived in Section 6. Three benchmark problems are solved in Section 7 to verify the method. Finally, we conclude our manuscript in Section 8.

2 Brief review of Maxwell equations

The general differential form of the Maxwell equations [Jin (2015)] are given by

with

E =electric field intensity (volts/meter),

D =electric flux density (coulombs/meter2),

H =magnetic field intensity (amperes/meter),

B =magnetic flux density (webers/meter2),

J =electric current density (amperes/meter2),

ρ=electric charge density (coulombs/meter2).

The divergence-free requirement in Eq. (1d) can be imposed for example with the penalty method [Rahman and Davies (1984)], vector finite elements [Whitney (2012); Nédélec(1980); Hano (1984)] or specially designed shape functions as presented in [Evans and Hughes (2013)]. The constitutive relations can be written as

where the constitutive parameters ?,μand σdenote, respectively, the permittivity(farady/meter), permeability (henrys/meter), and conductivity (siemens/meter) of the medium. For the time-harmonic fields with a single frequency, the time dependent parts of Maxwell’s equations can be written in simplified form as

where j is the imaginary unit and ωis the angular frequency. The vector wave equations forE can be obtained by eliminating Hfrom Eq. (3b) and considering the constitutive relations Eq. (2) to replaceH,

The boundary conditions for equations based on E are

3 Basic concepts in nonlocal operator method

Figure 1: (a) The electric field and notations. (b) Schematic diagram for support and dual-support, all circles above are support.Sx={x1, x2,x4,x6},Sx′={x1, x2,x3, x4}

Consider a domain as shown in Fig. 1(a), Let x be the spatial coordinates in the domain Ω;r:=x′-x is the spatial vector, the relative distance vector between x and x′;F:=F( x, t)and F′:=F( x′,t)are the electric field vectors for x and x′, respectively;Fr:=F′-F is the relative electric vector for vector r.

The vector wave equations can also be formulated by using onlyH. In this paper, we will employ the vector wave equations based on electric fields.

SupportSxis the domain where the nonlocal operator is defined, and any spatial point x′in support forms spatial vector r . Support Sxis usually presented by a spherical domain with a radius of hx. A point interacts with other points which fall inside the support of that point through nonlocal interactions.

In order to define the nonlocal operator, the shape tensor is defined as

which is symmetric. The prerequisites of shape tensor are that it shall be invertible, which can be satisfied usually when enough particles fall inside the support. Numerical example shows that the minimal number of neighbors in support is 2 and 3 for two-dimensional and three-dimensional problems, respectively.

Dual-supportis defined as the union of points whose supports include x , denoted by

Any point x′in Sx′forms a dual-vector r′(=-r). On the other hand,r′is the spatial vector formed inSx′. One example to illustrate the support and dual-support is shown in Fig. 1(b).

3.1 Nonlocal operators and definitions based on the support

In nonlocal operator method, key operators include the nonlocal operators for divergence,curl and gradient since they can be used to replace the local operators in the partial differential equations. We useto denote the nonlocal operator, while the local operator is ?. The nonlocal gradient of field F for point x in support Sxis defined as

with Fr=Fx′-Fx.

The nonlocal curl of field F for point x is defined as

The nonlocal divergence of field F for point x is defined as

The field value near a point x′can be approximated by Taylor series expansion by neglecting higher order terms as

Inserting Eq. (11) into the RHS of Eqs. (8), (9), (10) and integrating in support Sx, it can be shown that the nonlocal operators are identical to the local operators. For example,

When x′is close enough to x or when support xis small enough, the nonlocal operator can be considered as the linearization of the nonlinear field. The nonlocal operator converges to the local operator in the continuous limit. On the other hand, the nonlocal operator defined by integral form, still holds in the case where strong discontinuity exists and the local operator cannot be defined. The local operator can be viewed as a special case of the nonlocal operator.

3.2 Variation of nonlocal operators

The nonlocal operators defined above are in vector or tensor form. The variation of the nonlocal operators leads to a higher-order tensor form, which is not convenient for implementation. We need to express the high-order tensor into to vector or matrix form.Before we derive the variation of nonlocal operator, some notations to denote the variation and how the variations are related to the first- and second-order derivatives is to be discussed. Assuming a functional F (u, v), where u:=u(x),v:=v(x)are unknown functions in unknown vector[u, v], the first and second variation can be expressed as

It can be seen that the second variationis the double inner product of the Hessian matrix and the tensor formed by the variation of the unknowns, while the first variationis inner product of the gradient vector and the variation of the unknowns. The gradient vector and the Hessian matrix represent the residual vector and tangent stiffness matrix of the functional, respectively, with unknown functionsu, v being the independent variables,

The inner product or double inner product indicates that the location of an element in the residual or the tangent stiffness matrix corresponds to the location of the variation of the unknowns.

In this paper, we use a special variation

Obviously

The special first-order and second-order variation of a functional lead to the residual and tangent stiffness matrix directly. The traditional variation can be recovered by the inner product of the special variation and the variation of the unknown vector.

The number of dimensions of?δFxis infinite, and discretization is required.

After discretization of the domain by particles, the whole domain is represented by

where i is the global index of volume ΔVi,Nnode is the number of particles in ?.Particles inSiare represented by

where j1,..,jk,..,jniare the global indices of neighbors of particle i.

where denotes discretization,is all the variations of the unknowns in support Si,

i

where k is the index of particle jkin Ni. The process to obtainon nodal level is sometimes called the nodal assembly.

In the following section, we mainly discuss the special variation of the nonlocal operator and functional, while the actual variation can be recovered with ease.

where ΔVjkis the volume for particle jk. For the 3D case,is a 3× 3(ni+1)matrix, where niis the number of neighbors in Si,Niis given by Eq. (19). For each particle jk∈Nicalculating, we obtain

where k is the index of particle jk∈Ni. The minus sign denotes the reaction from the dual-support, which guarantees the regularity of the stiffness matrix in the absence of external constraints. The nodal assembly for the variation of the vector cross product can be finally obtained by

while the gradient of F×on {F0, F1, F2}is given by

The indices of R correspond to the locations in F×.

Similarly, the variation ofin the discrete form reads

where ΔVjkis the volumefor particle jk. In 3D,is a 9× 3(ni+1)matrix, where niis the number of neighbors in Si,Niis given by Eq. (19). For each particle in the neighbor list with, the terms in Rjkcan be added to theas

where k is the index of particle jk∈Ni. The sub-index ofcan be obtained by the way similar to Eq. (27).

4 Waveguide

A waveguide is a structure that guides waves, such as electromagnetic waves or sound,with minimal loss of energy by restricting expansions to one or two dimensions. The study of waveguide requires the eigenmode of the wave propagation which in turn requires the tangent stiffness matrix of the field. In this section, we derive the tangent stiffness matrix in the framework of variational principle using nonlocal operator method.The partial differential equation with boundary conditions for the waveguide problem can be expressed as

where Γ1is the electric boundary condition and Γ2is the magnetic boundary condition;?r(= ?/ ?0)and μr(=μ/μ0)denote the relative permittivity and relative permeability,respectively whilek0(=ω)is the wavenumber in free space,Z0(=)is the intrinsic impedance of free space and?0(=8.854× 10-12farad/meter) andμ0(=4π×10-7henry/meter) are the permittivity and permeability of the free space,respectively.

Consider the inner product of Eq. (30a) with arbitrary variationδEand integrate over the domain

Applying the same procedure to Eq. (30d) leads to

The sum of Eq. (31) and Eq. (32) is

Applying second vector Green's theorem in Eq. (34)

to Eq. (33) results in

Eq.35 is equivalent to the variation of the functional F (E),

The divergence-free condition is enforced by the penalty method, so the functional becomes

where p is the penalty parameter which is set to 1 in our examples as suggested in Rahman et al. [Rahman and Davies (1984)]. Finally, the eigenvalue problem of the waveguide problem reads

Eq. (38b) is the electric wall and Eq. (38c) is the magnetic wall, which is enforced for the sake of better accuracy and the elimination of some spurious solutions. For rectangular waveguide, the normal direction is parallel to a certain axis, for example n=(1,0,0),n× E =0?(Ey=0,Ez=0)Ey=0,Ez=0can be applied the same as Dirichlet boundary conditions.

F (E)on point x is

and its first variation is written as

Consider the first variation of all particles, and let, we have

The relation a?( b× c)=c?( a× b)=b?( c× a)is used in the third step. In the third and fourth steps of above derivation, the dual-support has been employed, i.e., in the third step, the termδFx′is the vector from x 's support, but is added to particle x′; since x′∈Sx,x belongs to the dual-supportof x′. In the fourth step, all the terms with δFxare collected from other particles whose supports contain xand therefore form the dual-support ofx . For any δFx, the first order variation δF(E)=0leads to the nonlocal form of the governing equation of the waveguide problem:

When the particle's volume ΔVx′→0, the continuous form is

Eq. (43) is the nonlocal governing equation of the waveguide on the electric field.For the eigenvalue problem, the stiffness matrix is required.

The special second variation of F (Ex)leads to the tangent stiffness matrix,

where

It should be noted that the latter termtakes the local value at the particle and can be obtained easily, while the nonlocal termcan be evaluated by Eq.(25). Assembling the stiffness matrix of all particles, one gets the global stiffness matrix and global “mass” matrix.

leading to the generalized eigenvalue problem

Note that the nodal integration of the above integrals results in hourglass modes which can be removed by introducing so-called hourglass energy, which will be addressed in the next section.

5 Hourglass energy functional

In order to remove the hourglass or zero-energy modes, a penalty term is added to achieve the linear completeness of the electric field, in which the penalty is proportional to the difference between the current value of a point and the value predicted by the field gradient.

The electric fieldin the neighborhood of a particle is required to be linear. Therefore, it has to be exactly described by the gradient of the electric field, and the hourglass modes are identified asthat part of the electric field, which is not described by the electric gradient. The difference between the current vector Erand the predicted vector by the field gradient (F:=in Eq. (8) is (F r-Er). We formulate the hourglass energy based on the difference in the support as follows: Letα=phg/(2μmK)be a coefficient for the hourglass energy, wheremK=tr(K),μis the magnetic coefficient,phgis the penalty which can be set to 1. Then, the functional for zero-energy mode is

The above definition of the hourglass energy is similar to the variance in probability theory and statistics. In the above derivation, we have used the relations,FTF: K=F:(F K), aTMb=M: a?b, A: B=tr(A BT), where the capital letter denotes the matrix and small letter the column vector. The purpose ofmKis to make the energy functional independent of the support since the shape tensorK is involved in FTF: K.

It should be noted that the zero-energy functional is valid in any dimensions and there is no limitation on the shape of the support.

Consider the variation of the zero-energy functional

The residual of the hourglass mode is required in the explicit time integration method. In this paper, we only discuss the implicit analysis.

The electric flux of the hourglass mode for one vector r is given by

Eq. (50) is an explicit formula for the hourglass flux. The term on δEis the hourglass term from its support, while the terms onδE′are the hourglass terms for the dual supportSx′′of point x′.

The second variation of the zero-energy functional is its stiffness matrix on one point.The global tangent stiffness matrix for hourglass energy functional can be assembled by

The above equations indicate, when the electric field is consistent with the field gradient,then the hourglass energy residual is zero.

Once the hourglass mode is eliminated, the residual of the hourglass mode is zero. The stiffness matrix of the hourglass mode overcomes the rank deficiency of the matrix system when nodal integration is used. The generalized eigenvalue problem becomes

where E is the eigenvector for all unknowns.

6 Nonlocal operator method for electromagnetic in the time domain

Consider a volume ?bounded by the surface S. The electromagnetic field generated by an electric current density Jxsatisfies the Maxwell equations. Eliminating the magnetic field with the aid of the constitutive relations, the curl-curl equation for the electric fieldE is obtained by Jin [Jin (2015)]:

In order to obtain the equivalent nonlocal form of Eq. (54), let us consider Eq. (40) from the previous section. Fromthe variational derivation of the waveguide problem,is equivalent to the functional

Based on Eq. (43), one can get the nonlocal form of

where

The Dirichlet boundary conditions are

where Pxis the specified electric wall on point x.

The Neumann boundary condition on S2can be written as

Finally, the central difference scheme can be used for the time integration yielding

7 Numerical examples

7.1 The Schr?dinger equation in 1D

In this section, we test the accuracy of the eigenvalue. The Schr?dinger equation written in adimensional units for a one-dimensional harmonic oscillator is

For simplicity, we use ω=1. The particles are distributed with constant/variable spacing Δxon the region [-10,10].

The exact wave functions and eigenvalues can be expressed as

where n is a non-negative integer.Hn(x)is the n-order Hermite polynomial.

The Schr?dinger equation in 1D is reformulated in variational form as

The tangent stiffness matrix is obtained as

We calculate the lowest eigenvalue and compare the numerical result with λ0=0.5. The convergence plot of the error with different weight function and discretizations is shown in Fig. 2. It can be seen that the convergence rate is around 2. The weight function and inhomogeneous discretization have limited effect on the convergence. Good agreement is observed between the numerical result and the exact solution.

Figure 2: Convergence of the lowest eigenvalue for a one-dimensional harmonic oscillator;1/r 3is the weight function; support is selected as h=nΔ x; dual-form with weight function 1/r3 uses an inhomogenous discretization in Fig. 3; the particle spacing in dual-form is selected as the minimal particle spacing in the discretization

The discretization of the dual-form is given in Fig. 3. The first three wave functions are given in Fig. 4.

Figure 3: The discretization of the dual form based on inhomogeneous discretization

Figure 4: The first three wave functions

7.2 Electrostatic field problems

When there is no electricity in the domain, Maxwell’s equations can be simplified into the Poisson equation with boundary conditions as

In the simulation, the boundary conditions are applied with the penalty method. The equivalent energy functional is

where α=1× 106is the penalty coefficient. The stiffness matrix can be obtained the similar way in Section 5.

In order to validate the accuracy of nonlocal operator formulation on the electrostatic field problem, we calculate the electrostatic field of a rectangular plate, which is a benchmark problem with the analytical solution given in Eq. (67). In this example, the potential on the upper and lower sides is 0, and the right side is, the horizontal electric field strength on the left side is 0, as shown in Fig. 5. In the simulation,the parameters are a=1m,b=1m,U0=1.0V.

Figure 5: Boundary condition of a rectangular plate

Figure 6: Contour plot of numerical solution of electric potential under mesh 50×50

Figure 7: Contour plot of analytical solution of electric potential

Figure 8: Convergence plot of L2 norm of electric potential

The support radius is selected as. The hourglass penalty of phg=1is used in all simulations. The plate is discretized with different mesh densities to test the convergence.The contour plot of the electric potential from the numerical solution and analytical solution are shown in Fig. 6 and Fig. 7, respectively. A satisfactory agreement is observed between Fig. 6 and Fig. 7. The L2 norm of the electric potential decreases with the refinement of the mesh with the convergence rate ofr=0.8709, as shown in Fig. 8.

7.3 Rectangular Waveguide problem

A hollow waveguide is a transmission line that looks like an empty metallic pipe. It supports the propagation of transverse electric (TE) and transverse magnetic (TM) modes,but not transverse electromagnetic (TEM) modes. There is an infinite number of modes that can propagate as long as the operating frequency is above the cutoff frequency of the mode. The notation TEmnand TMmnare commonly used to denote the type of the wave and its mode, where m and n are the mode number in the horizontal and vertical directions, respectively. The mode with the lowest cutoff frequency is called the fundamental mode or dominant mode. For a hollow rectangular waveguide, the dominant mode is TE10. The analytical solution for the E-field in the TE mode is expressed as

The electromagnetic analysis of a rectangular waveguide is well known [Pozar (2009)].Let us focus on the results used to verify our formulation, i.e.,

where fcmnis the cutoff frequency of mode mn ,kcmndenotes the wavenumber corresponding to mode mn while a and bare the width and height of the waveguide,respectively.

A section of a rectangular waveguide is modeled with the proposed nonlocal operator formulation and the first 3 modes are calculated and their field distributions analyzed.Since the background is set to a perfect electric conductor (PEC) material, we only need to model the vacuum inside the waveguide. The boundary conditions are “electric” in all directions, and the model is simulated using an eigenvalue solver in Matlab [Mathworks Guide (1998)]. In this model the first 3 modes are calculated. The dimensions of the waveguide are set to a=22.86mm,b=10.16mm and l=40mm; the boundaries in blue illustrate the electric walls and the red boundary is the magnetic wall, see Fig. 9. The domain of waveguide is discretized with two different particle spacings, as shown in Fig.10. The support is selected ash=2.2Δx and weight function w( r)=1/r2.

Figure 9: Section of a rectangular waveguide, where a=22.86 mm, b=10.16 mm and l=40 mm. Blue boundary denotes electric wall and red boundary is magnetic wall

Table 1: Comparison of fcmnbetween simulation and analytical results

The calculated frequencies are given in Tab. 1. The error in the frequency for Case 2 is less than 4%. Good agreements are obtained between the numerical results and theoretical results.The modes of the E-Field for two cases are shown in Fig. 11 and Fig. 12.

Figure 10: The discretizations for two cases

Figure 11: The TE modes for case 1

Figure 12: The TE modes for case 2

8 Conclusion

In this paper, we proposed a nonlocal operator formulation for electromagnetic problems employing variational principles. The formulation is implicit and provides the tangent stiffness matrix, which is needed for the solution of the eigenvalue problem. We presented a scheme for assembling the global stiffness matrix based on nonlocal operators. The nonlocal form of the electromagnetic vector wave equations based on the electric field is obtained by means of the variational principles. Three numerical examples, including the Schr?dinger equation in 1D, electrostatic field problem in 2D and waveguide in 3D are tested and show good agreement to available analytical solutions. In the future, we intend to solve also the transient problem and study problems involving strong discontinuities which are one of the key strength of nonlocal operator method.

Acknowledgement:The authors acknowledge the support of the German Research Foundation (DFG).