Scaled boundary isogeometric analysis for electrostatic problems

Engineering Analysis with Boundary Elements 85 (2017) 20–29

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

Scaled boundary isogeometric analysis for electrostatic problems Binghan Xue a,b,∗, Gao Lin a,b, Zhiqiang Hu a,b a b

State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian, Liaoning 116024, China Institute of Earthquake Engineering, Faculty of Infrastructure Engineering, Dalian University of Technology, Dalian, Liaoning 116024, China

a r t i c l e

i n f o

Keywords: Electrostatic problems SBFEM NURBS Scaled boundary isogeometric analysis Mortar method

a b s t r a c t The scaled boundary isogeometric analysis (SBIGA) is a novel semi-analytical technique, combing the advantages of the scaled boundary finite element method and the isogeometric analysis. In this paper, SBIGA is firstly exploited to solve electrostatic problems. According to the Laplace equation of electrostatic problems, the derivations and solutions of SBIGA equations for bounded domain and open domain problems are presented in details. A mortar method is employed to couple the solution on different subdomains, when the electrostatic problems with inhomogeneous media or complex boundaries which cannot be described by a single NURBS patch or cannot satisfy the scaling requirement in SBIGA. The mortar-based SBIGA can retain the flexibility of interface meshes compared with strong coupling methods. A condensation scheme is exploited to treat system equation in the analysis. Several numerical examples confirm the effectiveness, accuracy and convergence properties of SBIGA and the mortar-based SBIGA in solving electrostatic problems. © 2017 Elsevier Ltd. All rights reserved.

1. Introduction Electrostatic field analysis is one of the fundamental problems in computational electromagnetic, which disposes of the effects of electric charges at rest. For microelectromechanical systems (MEMS) and modern electron, an accurate electrostatic analysis is indispensable and essential. With the rapid development of the computer performance, numerical computational methods have been widely used in solving electrostatic problems. Numerical methods can provide remarkable results of the electrostatic problems, which may be difficult or even impossible to obtain by approaches of analytical calculations or experiments. Generally, such methods include finite difference method (FDM) [1], finite element method (FEM) [2,3], moment method [4,5], boundary element method (BEM) [6], meshless method [7,8] and scaled boundary finite element method (SBFEM) [9]. Among them, FDM which is based on replacing the differential quotient with the difference quotient approximately is the earliest used numerical method in solving electrostatic problems. In the framework of FDM, two dimensional domain and three dimensional domain are decomposed into equidistant or nonequidistant rectangular meshes and cubical meshes, respectively. Then, the difference equation of the system is obtained by replacing the differential equation respect to each node with the difference equation. But, the regular meshes cannot successfully simulate the electrostatic problem which has complicated material properties or complex geometries. FEM as the most commonly used numerical method is extensively exploited to solve the electrostatic problems. FEM has a remarkable ca-

∗

pability of simulating the electrostatic problems with complicated material properties and complex geometries. However, FEM is cumbersome in modelling problems with curved boundaries or open domain. Besides, performing mesh refinement and generating high quality meshes to fit the exact geometry is time consuming. BEM is also widely used in solving electrostatic problems. In BEM, only the boundary needs to be discretized, thus the spatial dimension is reduced by one and the data generation is easy. Furthermore, BEM is suitable to solve the open domain problems. However, a fundamental solution is generally complicated to express in a closed form for complex problems. SBFEM was proposed by Wolf and Song [10,11], which combines the merits of FEM and BEM, and has unique properties of its own. Like in FEM, the fundamental solution is unrequired. Like in BEM, only the boundary needs to be discretized, so the spatial dimension is reduced by one and the number of degrees of freedom is decreased. Establishing the scaled boundary coordinate system, that is, radial and circumferential direction, is one of the key steps in SBFEM. SBFEM is a semi-analytical method, since results are analytical in the radial direction and results are numerical in the circumferential direction. Moreover, SBFEM can handle the open domain problems by ranging the radial coordinate from one to infinite. However, like in FEM, the conventional Lagrange polynomials are chosen as the basis functions in SBFEM. Thus, the geometries with curved boundary cannot be exactly described and the continuity across element is C0 . In SBFEM, the scaling centre must be selected in a zone from which the whole boundary is visible, namely scaling requirement. When this condition

Corresponding author at: Faculty of Infrastructure Engineering, Dalian University of Technology, Dalian, Liaoning 116024, China. E-mail address: [email protected] (B. Xue).

https://doi.org/10.1016/j.enganabound.2017.09.012 Received 17 May 2017; Received in revised form 22 September 2017; Accepted 22 September 2017 0955-7997/© 2017 Elsevier Ltd. All rights reserved.

B. Xue et al.

Engineering Analysis with Boundary Elements 85 (2017) 20–29

is not satisfied, the total domain is generally subdivided in subdomains. The flexibility of the meshes is limited, since the meshes across the interfaces among the subdomains need to be point-wise matching. Isogeometric analysis (IGA) was proposed by Hughes et al. [12], in which the non-uniform rational B-splines (NURBS) are employed for geometries and field variables based on isoparametric concept. IGA has been widely used in thermal diffusion [13], contact problems [14], shape optimizing [15] and electromagnetics [16]. Compared with Lagrange polynomial-based numerical methods, IGA has some unique advantages. Firstly, CAD and CAE are integrated seamlessly by using the same NURBS basis functions. Secondly, the accuracy and robustness of IGA are superior, since the geometries with curved boundary can be exactly represented. Thirdly, any desired smoothness of NURBS basis function can be obtained by knot multiplicity. However, NURBS-based IGA has difficulty in treating the problems with complex topology, since a single NURBS patch cannot represent the complex topology. Thus, the complex domain is usually split into patches. Methods to couple the solution on different patches are needed. In order to retain the flexibility of the interface meshes, weak coupling methods are advantageous compared with strong point-wise couplings. Thus, the mortar method [17], which proposed for treating non-matching meshes in domain decomposing method, was employed as the coupling methods for different patches in IGA. The mortar-based IGA [18,19] can have an optical convergence and satisfy the inf-sup condition by choosing the appropriate Lagrange multipliers. More details about mortar methods can refer to literatures [20]. Inspiring by the notion of IGA, NURBS is also chosen as the basis functions in SBFEM. The NURBS-based SBFEM is named SBIGA, which has been successfully used in elastostatics problems [21] and fracture problems [22]. SBIGA fully combines the above-mentioned characteristics of SBFEM and IGA. Thus, SBIGA is firstly employed to solve the electrostatic problems in this paper. Moreover, a mortar-based SBIGA is proposed for treating the electrostatic problems with complex geometries, which cannot be described by a single NURBS patch or cannot satisfy the scaling requirement. In the mortar-based SBIGA, each patch can be represented and discretized by NURBS basis functions independently and the meshes across the interfaces are not forced to match point-wise any more. Besides, a condensation scheme which based on replace the degree of the freedom on the slave side with the degree of the freedom on the master side, is employed to simplify the system equation, then the condensed system equation can be solved with regular method in the analysis. This paper is structured as follows. In Section 2, an overview of Bsplines and NURBS basis functions is presented. In Section 3, the weak governing equation of electrostatic problems is derived by scaled boundary isogeometric analysis in detail. The mortar method and the condensation scheme are introduced in Section 4. In Section 5, some numerical examples are treated to verify the performance of SBIGA and the mortarbased SBIGA. The conclusions are presented in Section 6.

knot vector is called open knot vector which is the standard in CAD literature, and adopted in this paper. Univariate B-spline basis functions can be obtained by the Cox–De Boor formula [23]: { ⎧ 1, 𝑠 ∈ [𝑠𝑖 , 𝑠𝑖+1 ] ⎪𝑁𝑖,0 = 0, else (2) ⎨ 𝑠 −𝑠 ⎪𝑁𝑖,𝑝 = 𝑠−𝑠𝑖 𝑁𝑖,𝑝−1 (𝑠) + 𝑖+𝑝+1 𝑁𝑖+1,𝑝−1 (𝑠), 𝑝 ≥ 1 𝑠 − 𝑠 𝑠 − 𝑠 ⎩ 𝑖+𝑝 𝑖 𝑖+𝑝+1 𝑖+1 The important properties of B-spline basis functions are summarized below: ∑ a. Partition of unity: 𝑛𝑖=1 𝑁𝑖,𝑝 (𝑠) = 1, ∀𝑠 ∈ [𝑠1 , 𝑠𝑛+𝑝+1 ] b. Non-negativity: Ni,p (s) ≥ 0 c. Local support: Ni,p (s) = 0, ∀s∉[si ,si + p + 1 ) ∑ d. Linear independence: 𝑛𝑖=1 𝛼𝑖 𝑁𝑖,𝑝 (𝑠) = 0 ⇔ 𝛼𝑖 = 0 e. Flexible continuity: a knot with multiplicity k means the smoothness of the basis function is Cp-k at that knot. If a knot with multiplicity p, the basis function is C0 continuity and has interpolatory at that knot. 2.2. NURBS basis functions NURBS basis functions are defined by the B-splines and corresponding weights. Univariate NURBS basis functions can be given by: 𝑅𝑝𝑖 (𝑠) = ∑𝑛

𝜔𝑖 𝑁𝑖,𝑝 (𝑠)

𝑗=1

𝜔𝑗 𝑁𝑗,𝑝 (𝑠)

(3)

where 𝜔i are the corresponding weights of B-spline basis functions Ni,p . NURBS basis functions inherit the feature of the B-spline basis functions, such as partition of unity, local support, non-negativity, linear independence and flexible continuity. Comparing with B-spline basis functions, NURBS basis functions can exact construct conic sections such as circles and ellipses. A NURBS curve is constructed as a linear combination of NURBS basis functions: 𝑛 ∑ 𝐶 (𝑠 ) = 𝑅𝑝𝑖 (𝑠)𝐏𝑖 (4) 𝑖=1

Fig. 1 shows quadratic NURBS basis functions respect to the knot vector s = {0,0,0,1,2,2,3,4,4,4}. The corresponding NURBS curve is shown in Fig. 2. 3. The scaled boundary isogeometric analysis The scaled boundary isogeometric analysis is a new method by combining the isogeometric analysis and scaled boundary finite element method. In this section, the weak governing equation of electrostatic

2. B-splines and NURBS basis functions The basic concepts of B-splines and NURBS basis functions, which will be used in SBIGA, are briefly introduced in this section. For more details, we can refer to the literature [23,24]. 2.1. B-splines A B-spline basis is denoted by a non-decreasing sequence of real numbers: { } 𝐬 = 𝑠1 , 𝑠2 , … , 𝑠𝑛+𝑝+1 (1) where s is called a knot vector, each si is a knot, n is the number of basis functions and p is the degree of basis functions. Moreover, [s1 , sn+p+ 1 ] defines a single B-spline patch and the interval [si , si +1 ] is the ith knot span. If the first and last entries in s have a multiplicity of p + 1, this

Fig. 1. Quadratic NURBS basis functions.

21

B. Xue et al.

Engineering Analysis with Boundary Elements 85 (2017) 20–29

The gradient operators of these two coordinate systems are related as: ⎧ 𝜕 ⎫ ⎧ 𝜕 ⎫ ⎪ 𝜕 𝑥̂ ⎪ −1 ⎪ 𝜕𝜉 ⎪ = 𝐉 (9) ⎨ 𝜕 ⎬ ⎨1 𝜕 ⎬ ⎪ ⎪ ⎪ ⎪ ⎩ 𝜕 𝑦̂ ⎭ ⎩ 𝜉 𝜕𝜂 ⎭ where, J is the boundary Jacobian matrix which is independence of 𝜉: [ ] 𝑥 𝑦 𝐉(𝜂) = (10) 𝑥,𝜂 𝑦,𝜂 Writing the derivative of J as the following form: [ ] [ ] 𝑗 𝑗12 𝐉−1 (𝜂) = 11 = 𝒃1 (𝜂) 𝒃2 (𝜂) 𝑗21 𝑗22

Fig. 2. Quadratic NURBS curve.

(11)

Then, the Laplace operator can be written as: ∇ = 𝐛1 (𝜂)

𝜕 1 𝜕 + 𝐛 (𝜂) 𝜕𝜉 𝜉 2 𝜕𝜂

(12)

The infinitesimal area is expressed as: 𝑑Ω = 𝜉|𝐉(𝜂)|𝑑𝜉𝑑𝜂

(13)

3.3. Boundary isogeometric discretization In the scaled boundary coordinate system, the boundary of the domain Ω includes and outer boundary Γb and side boundary Γs . In the framework of SBIGA, only the outer boundary is represented and isogeometrically discretized by the NURBS basis functions introduced in Section 2.2. Unlike numerical methods based on the Lagrange basis functions, deviation is inexistence between the geometry model and the computing model. An arbitrary point (x(𝜂), y(𝜂)) on the outer boundary is represented by NURBS and expressed as:

Fig. 3. The coordinate definition and isogeometric discretization of SBIGA.

problems is derived by using scaled boundary isogeometric analysis in detail. 3.1. Basic equations in electrostatic field

𝑥(𝜂) = 𝐑(𝜂)𝐱 𝑦(𝜂) = 𝐑(𝜂)𝐲

The basic equations in electrostatic field problems can be written as Laplace equation: ∇2 𝜀𝜙 = 0,

in 𝛀

where the circumstance coordinate 𝜂 belongs to the parameter space spanned by the knot vector s, (x, y) denotes the Cartesian coordinate of control points corresponding to NURBS basis function R(𝜂). Substituting Eq. (14) into Eq. (8), then the coordinate of an arbitrary point in the domain can be rewritten as:

(5)

with essential boundary conditions ̄ 𝜙 = 𝜙,

on 𝚪𝐷

(6)

and natural boundary conditions 𝜙,𝑛 = 𝑞̄,

on 𝚪𝑁

(14)

𝑥̂ = 𝜉(𝐑(𝜂)𝐱 − 𝑥̂ 0 ) + 𝑥̂ 0 𝑦̂ = 𝜉(𝐑(𝜂)𝐲 − 𝑦̂0 ) + 𝑦̂0

(7)

(15)

where 𝜑 denotes potential, 𝜀 denotes permittivity. 𝜙̄ and 𝑞̄ denote the given potential and electric intensity normal to the boundary.

Based on the isoparametric concept, the potential field is also represented by the NURBS basis functions. Then, the potential of an arbitrary point in the domain can be written as:

3.2. Coordinate system transformation

𝜙(𝜉, 𝜂) = 𝐑(𝜂)𝝓(𝜉)

In order to apply SBIGA to solve 2D electrostatic problems, as shown in Fig. 3, two coordinate systems are established in SBIGA. That is, the Cartesian coordinate system (x, y) and the scaled boundary coordinate system (𝜉, 𝜂). In the scaled boundary coordinate system, scaling centre O needs to be chosen in a zone which satisfies the scaling requirement. 𝜉 is the radial coordinate. For bounded domain, 𝜉 is ranging from 0 at the scaling centre to 1 on the outer boundary Γb . For unbounded domain, 𝜉 is ranging from 1 on the outer boundary Γb to infinity. 𝜂 is the circumferential coordinate. The mapping relationship between the Cartesian coordinate system and the scaled boundary coordinate system is formulated as:

where 𝝋(𝜉) denotes the potential at control points which is analytical in radial direction. 𝝋(𝜉 = 1) represents potential on the outer boundary. 3.4. Derivation of the SBIGA governing equations The equivalent variational problem of the electrostatic problem Eqs. (5)–(7) is defined as: 𝛿Π =

∫Ω

𝜀(∇𝛿𝜙)𝑇 ∇𝜙𝑑Ω −

∫Γ𝑏𝑁

𝛿𝜙(𝜂)𝑇 𝑞̄𝑏 (𝜂)𝑑Γ−

∫Γ𝑠𝑁

𝛿𝜙(𝜂)𝑇 𝑞̄𝑠 (𝜉)𝑑𝜉

(17)

where, 𝑞̄𝑏 and 𝑞̄𝑠 are the normal electric intensity respect to outer boundary Γb and side boundary Γs , respectively. Substituting Eqs. (12), (13) and (16) into Eq. (17), 𝜉1 [ 1 𝛿𝝓T (𝐄0 𝝓,𝜉 + (𝐄1 )T 𝝓 − 𝐏)|𝜉=1 − 𝛿𝝓(𝜉)T (𝐄0 𝜉 2 𝝓(𝜉),𝜉𝜉 ∫𝜉0 𝜉 ] + (𝐄0 + (𝐄1 )T − 𝐄1 )𝜉𝝓(𝜉),𝜉 − 𝐄2 𝝓(𝜉) + 𝐅(𝜉)) 𝑑𝜉 = 0 (18)

𝑥̂ (𝜉, 𝜂) = 𝜉(𝑥(𝜂) − 𝑥̂ 0 )) + 𝑥̂ 0 𝑦̂(𝜉, 𝜂) = 𝜉(𝑦(𝜂) − 𝑦̂0 )) + 𝑦̂0

(16)

(8)

where (𝑥̂ , 𝑦̂) represents the Cartesian coordinate of an arbitrary point in the physical domain. (𝑥̂ 0 , 𝑦̂0 ) signifies the scaling centre O. (x(𝜂), y(𝜂)) denotes the Cartesian coordinate of an arbitrary point on the outer boundary.

with 22

B. Xue et al.

𝐄0 = 𝐄1 = 𝐄2 = 𝐏=

∫Γ𝑏 ∫Γ𝑏 ∫Γ𝑏

Engineering Analysis with Boundary Elements 85 (2017) 20–29

𝜀𝐁1 (𝜂)T 𝐁1 (𝜂)|𝐉|𝑑𝜂 𝜀𝐁2 (𝜂)T 𝐁1 (𝜂)|𝐉|𝑑𝜂 𝜀𝐁2 (𝜂)T 𝐁2 (𝜂)|𝐉|𝑑𝜂

∫Γ𝑏𝑁

𝐅(𝜉) = 𝜉

𝐑(𝜂)T 𝑞̄𝑏 (𝜂)𝑑Γ

∫Γ𝑠𝑁

𝑞̄𝑠 (𝜉)𝑑Γ

(19)

where, E0 , E1 and E2 are the constant coefficient matrix. P is the equivalent electric charge at control points. F is the electric charge due to electric intensity respect to side boundary. Considering the arbitrariness of 𝛿𝝋 and 𝛿𝝋(𝜉), then Eq. (18) can be rewritten as: 𝐄0 𝜉 2 𝝓(𝜉),𝜉𝜉 + (𝐄0 + (𝐄1 )T − 𝐄1 )𝜉𝝓(𝜉),𝜉 − 𝐄2 𝝓(𝜉) + 𝐅(𝜉) = 0

Fig. 4. The sketch of domain decomposition.

(20) Eliminating the coefficient c1 in Eq. (29) and considering Eq. (24) lead to:

with 𝐄0 𝝓,𝜉 + (𝐄1 )T 𝝓 = 𝐏

(21)

𝐊𝝓 − 𝐏 = 0

Eq. (20) is the weak governing equation of the system and can be solved analytically. Eq. (21) is the corresponding boundary conditions on the boundary 𝜉 = 1.

(30)

with 𝐊 = 𝚽21 (𝚽11 )−1

(31)

For an unbounded domain, the potential at 𝜉 = ∞ must be finite, thus the coefficient c1 = 0. Then Eq. (28) can be rewritten as: { 𝝓(𝜉) = 𝚽12 𝜉 −𝛃𝑝 𝐜2 (32) 𝐪(𝜉) = 𝚽22 𝜉 −𝛃𝑝 𝐜2

3.5. Solution scheme for the SBIGA governing equations When no normal electric intensity acts on the side boundary, that is F = 0, Eq. (20) degenerates to homogeneous second-order differential equation, which can be transformed to two first-order differential equations with two unknowns. Introducing a new variable: { } 𝝓(𝜉) 𝐗(𝜉) = (22) 𝐪(𝜉)

Eliminating the coefficient c2 in Eq. (32) and considering Eq. (24) lead to: 𝐊𝝓 − 𝐏 = 0

(33)

with 𝐊 = −𝚽22 (𝚽12 )−1

where 𝐪(𝜉) = 𝐄0 𝜉𝝓(𝜉),𝜉 + (𝐄1 )T 𝝓(𝜉)

Eqs. (30) and (33) are algebraic equations on the boundary 𝜉 = 1 for the bounded domain and the unbounded domain, respectively. The boundary potential can be calculated by Eq. (30) or (33). Then, the potential inside the domain is obtained by the first line in Eq. (28).

(23)

On the outer boundary 𝜉 = 1, Eq. (23) can be written as: 𝐪(𝜉 = 1) = 𝐄0 𝝓(𝜉 = 1),𝜉 + (𝐄1 )T 𝝓(𝜉 = 1) = 𝐏

(34)

(24)

Substituting Eqs. (22) and (23) into Eq. (20) without the nonhomogeneous term leads to:

4. The mortar method

𝜉𝐗(𝜉),𝜉 = −𝐙𝐗(𝜉)

For electrostatic problems, the solving process of SBIGA is introduced in Section 3. When the geometry of the computational domain is complex, which cannot be described by a single NURBS patch or cannot satisfy the scaling requirement, the total domain must be decomposed into subdomains. Each subdomain can be described by a single NURBS patch and satisfy the scaling requirement. Thus, methods to couple the numerical solutions on different patches are required. In the framework of SBIGA, a mortar method is introduced as the coupling method in this section.

(25)

with a Hamiltonian matrix: [ ] (𝐄0 )−1 (𝐄1 )T −(𝐄0 )−1 𝐙= 𝐄1 (𝐄0 )−1 (𝐄1 )T − 𝐄2 −𝐄1 (𝐄0 )−1

(26)

Eq. (25) is usually solved by eigenvalue method [11]. The eigenvectors 𝚽 and eigenvalue 𝚲 of the Hamiltonian matrix can be formulated as: [ ][ ] 𝚽11 𝚽12 𝛃𝑛 𝐙𝚽 = 𝚽𝚲 = (27) 𝚽21 𝚽22 𝛃𝑝 {

4.1. Domain decomposition As shown in Fig. 4, the domain Ω is decomposed into K nonoverlapping subdomains Ωk . Each subdomain is represented by a single NURBS patch and satisfies the scaling requirement. Topology relationships among these subdomains are formulated as:

Substituting Eq. (27) into Eq. (25) leads to: 𝝓(𝜉) = 𝚽11 𝜉 −𝛃𝑛 𝐜1 + 𝚽12 𝜉

−𝛃𝑝

𝐜2

𝐪(𝜉) = 𝚽21 𝜉 −𝛃𝑛 𝐜1 + 𝚽22 𝜉 −𝛃𝑝 𝐜2

(28)

For a bounded domain, the potential at 𝜉 = 0 must be finite, thus the coefficient c2 = 0. Then Eq. (28) can be rewritten as: { 𝝓(𝜉) = 𝚽11 𝜉 −𝛃𝑛 𝐜1 (29) 𝐪(𝜉) = 𝚽21 𝜉 −𝛃𝑛 𝐜1

̄ = Ω

𝐾 ⋃ 𝑘=1

̄𝑘 Ω

Ω𝑘 ∩ Ω𝑙 = ∅ 𝜕 Ω𝑘 ∩ 𝜕 Ω𝑙 = Γ𝑘𝑙 , 23

𝑘≠𝑙

(35)

B. Xue et al.

Engineering Analysis with Boundary Elements 85 (2017) 20–29

̄ 𝑘, where, Ωk and 𝜕Ωk are the interior and boundary of subdomain Ω ̄ 𝑘 and subdomain respectively. Γkl is the interface between subdomain Ω ̄ 𝑙. Ω The continuity condition across the interface can be satisfied by the mortar weak coupling method: ∑𝑁 𝑖=1

∫Γ𝑖

[𝜙]𝑖 𝜆𝑑Γ = 0

4.3. The mortar integration scheme An element based integral scheme [26,27] is exploited to calculate the mortar integral equation (40). In the element based integral scheme, the quadrature point is chosen as the Gauss integral points for each element over the slave side. This scheme is rather efficient, because no information from the master side is included when determining the locations of the Gauss points. According to the Gauss quadrature rule, one element of the matrix GsΓ can be calculated as:

(36)

where N is number of interfaces, [𝜑]i is the jump across the interface Γi . 𝜆 is the Lagrange multiplier and gives an approximation of the normal flux across the interface. For conciseness, the number of interfaces N is set to be one, that is, the total domain is decomposed into two subdomains. One of the subdomains with refined meshes is chosen as the slave (non-mortar) side Ωs and one with coarse meshes is chosen as the master (mortar) side Ωm in generally. Each subdomain is discretized by the NURBS basis functions, independently. The jump [𝜑] is calculated as: [𝜙]=𝐑𝑠Γ 𝝓𝑠Γ − 𝐑𝑚Γ 𝝓𝑚Γ

𝐺𝑠𝑖𝑘Γ =

where subscripts sΓ and mΓ denote interface parameters and variables respect to the slave side and the master side, respectively. The Lagrange multiply is typically defined over the slave side with an equivalent discretization:

𝑖𝑘 𝐺𝑚 Γ =

̂ is the basis functions in the Lagrange multiplier space, 𝝀 is where 𝐑 Lagrange multipliers on control points. Substituting Eqs. (37) and (38) into Eq. (36) (39)

with 𝐆𝑠Γ = 𝐆𝑚Γ =

∫Γ ∫Γ

̂ 𝜆 𝐑𝑠Γ 𝑑Γ 𝐑 ̂ 𝜆 𝐑𝑚Γ 𝑑Γ 𝐑

𝑙=𝑙

( ) ( ) 𝑅̂ 𝑖𝜆 𝜂𝑠𝑙 Γ 𝑅𝑘𝑠Γ 𝜂𝑠𝑙 Γ ||𝐽𝑙 ||𝑊𝑙

(43)

𝑁int

∑ 𝑙=𝑙

( ) ( ) 𝑅̂ 𝑖𝜆 𝜂𝑠𝑙 Γ 𝑅𝑘𝑚Γ 𝜂𝑚𝑙 Γ ||𝐽𝑙 ||𝑊𝑙

(44)

Noteworthy, 𝑅𝑘𝑚Γ (𝜂𝑚𝑙 Γ ) and 𝑅̂ 𝑖𝜆 (𝜂𝑠𝑙 Γ ) are defined on the master side and the salve side, respectively. The integral point 𝜂𝑚𝑙 Γ on the master side must respect to the integral point 𝜂𝑠𝑙 Γ on the slave side. Thus, calculating the matrix GmΓ is more complex than that of the matrix GsΓ . A fundamental operation of NURBS, namely point inversion algorithm [23] is employed to determine 𝜂𝑚𝑙 Γ which corresponds to 𝜂𝑠𝑙 Γ . The point inversion algorithm which belongs to iterative methods can be expressed as: ( ) ( ( ) ( )) 𝐶,𝜂 𝜂𝑚𝑙,𝑖Γ ⋅ 𝐶 𝜂𝑚𝑙,𝑖Γ − 𝐶 𝜂𝑠𝑙 Γ 𝑙,𝑖+1 𝑙,𝑖 𝜂𝑚Γ = 𝜂𝑚Γ − (45) ( ) ( ( ) ( )) | ( )|2 𝐶,𝜂𝜂 𝜂𝑚𝑙,𝑖Γ ⋅ 𝐶 𝜂𝑚𝑙,𝑖Γ − 𝐶 𝜂𝑠𝑙 Γ + ||𝐶,𝜂 𝜂𝑚𝑙,𝑖Γ || | | where C, 𝜂 and C, 𝜂𝜂 mean the first derivative and the second derivative of the interface on the master side. The point inversion algorithm is performed for each integral point on the master side.

(38)

𝐆 𝑠 Γ 𝝓𝑠 Γ = 𝐆 𝑚 Γ 𝝓𝑚 Γ

∑

where Nint is the number of integral points. Jl and Wl are the Jacobian and the integral weight on the integral point l, respectively. 𝜂𝑠𝑙 Γ is the integral point on the salve side. 𝑅𝑘𝑠Γ and 𝑅̂ 𝑖𝜆 are determined by 𝜂𝑠𝑙 Γ and denote the NURBS basis function and the Lagrange multiplier basis function, respectively. The matrix GsΓ can be obtained by looping the integration equation (43) on each element over the salve side. Similarly, one element of the matrix GmΓ can also be calculated as:

(37)

̂ 𝜆𝝀 𝜆=𝐑

𝑁int

(40)

Eq. (39) is continuity conditions across the interface, which will be used to implement the condensation scheme. Eq. (40) is the mortar integral equation.

4.4. The condensation scheme 4.2. The Lagrange multiplier space By using the scaled boundary isogeometric analysis, the algebraic equation of subdomains Ωs and Ωm can be combined as: [ ][ ] [ ] 𝐊𝑠 𝝓𝑠 𝐏𝑠 = (46) 𝐊 𝑚 𝝓𝑚 𝐏𝑚

In order to guarantee the mortar method of optimal order and to be well posed, two abstract requirements should be satisfied [25]. These requirements include an appropriate approximate order of the Lagrange multiplier space and a uniform inf-sup stability between the primal space and the Lagrange multiplier space. In this paper, we employ the NURBS basis function with a cross point modification to construct the Lagrange multiplier space. This construction method (p/p pairing) [18] has been proven to have an optical convergence order p and sat̂ in the Lagrange isfy the uniform inf-sup stability. The basis function 𝐑 multiplier space is defined as: 𝑝 𝑝 𝑖 ∈ {2, … , 𝑝 + 1} ⎧𝑅𝑖 (𝜂) + 𝛼𝑖 𝑅1 (𝜂), ⎪ 𝑝 𝑝 ̂ 𝑅𝑖 (𝜂) = ⎨𝑅𝑖 (𝜂), 𝑖 ∈ {𝑝 + 2, … , 𝑛 − 𝑝 − 1} ⎪ 𝑝 𝑝 ⎩𝑅𝑖 (𝜂) + 𝛽𝑖 𝑅𝑛 (𝜂), 𝑖 ∈ {𝑛 − 𝑝, … , 𝑛 − 1}

Rewriting the potential 𝝋s and 𝝋m as: [

𝛽𝑖 = −𝑅𝑝𝑖 (𝑝) (𝜂)∕𝑅𝑝𝑛(𝑝) (𝜂),

( ) 𝜂 ∈ 0, 𝜂2 ( ) 𝜂 ∈ 𝜂𝐸−1 , 1

(47)

where subscripts sΓ and mΓ denote interface parameters and variables respect to the slave side and the master side, respectively. Subscripts si and mi denote the rest parameters and variables respect to the slave side and the master side, respectively. Substituting Eq. (39) into Eq. (47) leads to:

(41)

with 𝛼𝑖 = −𝑅𝑝𝑖 (𝑝) (𝜂)∕𝑅𝑝1(𝑝) (𝜂),

⎡𝝓 ⎤ ] ⎢ 𝑠Γ ⎥ 𝝓𝑠 𝝓 = ⎢ 𝑠𝑖 ⎥ ⎢𝝓𝑚Γ ⎥ 𝝓𝑚 ⎢𝝓 ⎥ ⎣ 𝑚𝑖 ⎦

[ (42)

̂ In the Lagrange multiplier space, the degree of the basis functions 𝐑 is reduced by one on the element corresponding to the cross point, and the inter-element continuity is retained.

⎡ 𝝓 ⎤ ⎡0 ] ⎢ 𝑠Γ ⎥ ⎢ 𝝓𝑠 𝝓 𝐈 = ⎢ 𝑠𝑖 ⎥ = ⎢ ⎢𝝓𝑚Γ ⎥ ⎢0 𝝓𝑚 ⎢ 𝝓 ⎥ ⎢0 ⎣ 𝑚𝑖 ⎦ ⎣

with ( )−1 𝐌 = 𝐆𝑠Γ 𝐆𝑚Γ 24

𝐌 0 𝐈 0

0⎤ ⎡ 𝝓𝑠𝑖 ⎤ ⎥⎡ 𝝓 ⎤ 0⎥⎢ 𝑠𝑖 ⎥ ̂ ⎢𝝓 ⎥ 𝝓𝑚 Γ = 𝐌 ⎥ ⎢ 𝑚Γ ⎥ 0⎥⎢ ⎣ 𝝓𝑚𝑖 ⎦ ⎣ 𝝓𝑚𝑖 ⎦ ⎥ 𝐈⎦

(48)

(49)

B. Xue et al.

Engineering Analysis with Boundary Elements 85 (2017) 20–29

Fig. 5. The model of the first example.

Noteworthy, GsΓ is an invertible matrix by using the constructed Lagrange multiplier space. Substituting Eq. (48) into Eq. (46): ⎡𝐊𝑠ΓΓ ⎢ ⎢ 𝐊𝑠𝑖Γ ⎢ ⎢ ⎣

𝐊𝑠Γ𝑖 𝐊𝑠𝑖𝑖

⎤ ⎡ 𝐏𝑠Γ ⎤ ⎥ ⎡ 𝝓𝑠𝑖 ⎤ ⎢ ⎥ ⎥𝐌 ̂ ⎢𝝓 ⎥ = ⎢ 𝐏𝑠𝑖 ⎥ 𝑚Γ ⎥ ⎢𝐏 ⎥ 𝐊𝑚Γ𝑖 ⎥ ⎢ ⎣ 𝝓𝑚𝑖 ⎦ ⎢ 𝑚Γ ⎥ 𝐊𝑚𝑖𝑖 ⎥⎦ ⎣ 𝐏𝑚𝑖 ⎦

𝐊𝑚ΓΓ 𝐊𝑚𝑖Γ

(50)

Fig. 6. Comparison of the exactness of the boundary representation.

̂ Then, left multiply both sides of Eq. (50) by 𝐌 ⎡ 𝐊𝑠𝑖𝑖 ⎢ 𝑇 ⎢𝐌 𝐊𝑠Γ𝑖 ⎢ ⎣ 0

𝐊𝑠𝑖Γ 𝐌 𝐌𝑇 𝐊

𝑠ΓΓ 𝐌

+ 𝐊𝑚ΓΓ

𝐊𝑚𝑖Γ

0 ⎤⎡ 𝝓 ⎤ ⎡ 𝐏𝑠𝑖 ⎤ ⎥ 𝑠𝑖 𝐊𝑚Γ𝑖 ⎥⎢𝝓𝑚Γ ⎥ = ⎢𝐌𝑇 𝐏𝑠Γ + 𝐏𝑚Γ ⎥ ⎢ ⎥ ⎢ ⎥ ⎥ 𝐏𝑚𝑖 ⎦ 𝐊 ⎦⎣ 𝝓𝑚𝑖 ⎦ ⎣

Table 1 Information of different meshes.

(51)

Mesh

Number of elements

𝑚𝑖𝑖

Eq. (51) is the condensate system equation of the total domain Ω. The coefficient of Eq. (51) is a sparse symmetric and positive definite matrix. Then, the boundary potential of each subdomain can be obtained by solving Eq. (51).

1 2 3 4 5

8 16 32 64 128

DOFs SBIGA

SBFEM

12 20 36 68 132

16 32 64 128 256

5. Numerical examples In this section, we will investigate the performance of the scaled boundary isogeometric analysis and the mortar-based scaled boundary isogeometric analysis in solving electrostatic problems. To check the convergence of SBIGA and mortar-based SBIGA, the root-mean-square (RMS) error is defined as: √ ∑( )2 ∑ ( )2 RMS = (52) 𝜙𝑖_exact − 𝜙𝑖_calc ∕ 𝜙𝑖_exact 𝑖

𝑖

where, 𝜙𝑖_exact is analytical solution, 𝜙𝑖_calc is the numerical solution. 5.1. Infinitely long dielectric cylinder An infinitely long dielectric cylinder with natural boundary condition is presented to verify the performance of SBIGA to analyse the problems with curved boundary. For conciseness, one section of the cylinder as shown in Fig. 5 is considered in the analysis. The section can be described by a single NURBS patch and satisfy the scaling requirement, thus there is no need to use the mortar method. In the analysis, the centre of the section is chosen as the scaling centre O and SBFEM is also employed as a comparison. In the analysis, the domain boundary is discretized by the quadratic Lagrange basis function in SBFEM and the quadratic NURBS basis function in SBIGA, respectively. The meshes consist of eight elements with the same density was obtained by different discretization scheme. Fig. 6 presents the distance between the point on the boundary and the scaling centre O. The exact distance is the radius 1. It is observed that the NURBS basis function can exactly represent the boundary pointwisely, while the Lagrange basis function describes the boundary approximately except the nodal location. The analytical solution of this example is expressed as: 𝜙(𝑟, 𝜃) = 𝑟0 cos 𝜃

1

Fig. 7. RMS error of the potential on the periphery.

and SBFEM. The RMS errors of the potential on the periphery for SBIGA and SBFEM are presented in Fig. 7. As shown in Fig. 7, SBIGA has a more rapid convergence rate than SBFEM with the increasing of element number. The comparison of RMS error versus CPU time is presented in Fig. 8. The CPU time is recorded in seconds on a Windows 7 desktop computer (CPU: Intel core i7-3770, 3.40 GHz). As shown in Fig. 8, SBIGA consumes less CPU time than SBFEM when achieves the same level accuracy. Thus, the accuracy and efficiency of SBIGA are better than that of SBFEM. 5.2. Infinitely long square metal trough with inhomogeneous media

(53)

In order to verify the effectiveness of the mortar-based SBIGA in solving electrostatic problems, an infinitely long square metal trough problem with inhomogeneous media which has an analytical solution

Five meshes with increasing element number as shown in Table 1 are exploited to compare the convergence and accuracy between SBIGA 25

B. Xue et al.

Engineering Analysis with Boundary Elements 85 (2017) 20–29

Fig. 10. RMS errors of the potential on the total domain (P = 2). Fig. 8. RMS error versus CPU time.

Ω2 . O1 and O2 are chosen as the scaling centre in subdomains Ω1 and Ω2 , respectively. In the mortar-based SBIGA, each subdomain is represented by a single NURBS patch. A foundational operation of NURBS basis function, namely h-refinement algorithm is employed to discretize each subdomain independently. Thus, subdomains Ω1 and Ω2 have different element density and the meshes across the interface Γ are nonmatching. One of the subdomains with much more refined meshes is chosen as the slave side and another one is chosen as the master side. Four different meshes with increasing element number are obtained by using the h-refinement algorithm. A strong coupling method is also employed as a comparison. In the strong coupling method, the meshes across the interface must be matching pointwisely and four different meshes with increasing element number are also used. The analytical solution of this problem is expressed as: { [ ] sinh (𝑘𝑥) 𝑎1 sinh (𝑘𝑦) + 𝑏1 cosh (𝑘𝑦) 𝑐≤𝑦≤𝑏 𝜙(𝑥, 𝑦) = (54) 𝑐1 sin (𝑘𝑥) sinh (𝑘𝑦) 0≤𝑦≤𝑐

Fig. 9. Sketch for the rectangular metal trough.

is presented. The geometry and the boundary conditions of the square metal trough are presented in Fig. 9. In the analysis: a = 1, b = 1, c = 0.5, U0 = 100, 𝜀2 = 1.25𝜀1 . The domain Ω is decomposed into two subdomains Ω1 and Ω2 , then Γ is the interface between subdomains Ω1 and

with 𝑘 = 𝜋∕𝑎 [ ( ) ] 𝑑1 = sinh (𝑘𝑏) 𝜀1 tanh (𝑘𝑐 ) − 𝜀2 coth (𝑘𝑐 ) + 𝜀2 − 𝜀1 coth (𝑘𝑏)

Fig. 11. Distribution of the electric intensity in X-direction.

26

B. Xue et al.

Engineering Analysis with Boundary Elements 85 (2017) 20–29

Fig. 12. Distribution of the electric intensity in Y-direction.

Fig. 13. Model of the U-shaped metal trough.

[ ] 𝑎1 = 𝑈0 𝜀1 tanh (𝑘𝑐 ) − 𝜀2 coth (𝑘𝑐 ) ∕𝑑1 ( ) 𝑏1 = 𝑈0 𝜀2 − 𝜀1 ∕𝑑1 [ ( ) ] 𝑐1 = 𝑈0 𝜀1 tanh (𝑘𝑐 ) − 𝜀2 coth (𝑘𝑐 ) + 𝜀2 − 𝜀1 coth (𝑘𝑏) ∕𝑑1

sis functions, independently. Subdomain Ω1 is subdivided into 128 elements and chosen as the slave side, and subdomain Ω2 is subdivided into 64 elements and chosen as the master side. The meshes across the interface Γ are non-matching. Numerical analysis using ANSYS is employed to provide the reference results. In ANSYS, the U-shaped metal trough is subdivided into 2800 elements. As shown in Fig. 14, the distribution of potential obtained by the mortar-based SBIGA is nearly identical compared with that of ANSYS. Thus, the mortar-based SBIGA is powerful method to solve the electrostatic problems with complex boundary.

(55)

As shown in Fig. 10, the RMS errors of the mortar-based SBIGA with non-matching meshes and SBIGA with matching meshes are nearly equal and have an approximate slope. As shown in Figs. 11 and 12, the distribution of the electric intensity in X and Y directions calculated by the mortar-based SBIGA is identical with that of the analytical solution. Thus, the mortar-based SBIGA is effective and accurate in solving electrostatic problem with inhomogeneous media. 5.3. Infinitely long U-shaped metal trough

5.4. Open domain problem

As shown in Fig. 13(a), the infinitely long U-shaped metal trough cannot satisfy the scaling requirement, when represented by a single NURBS patch. Thus, as shown in Fig. 13(b), the total domain Ω needs to be decomposed into subdomains Ω1 and Ω2 .O1 and O2 are chosen as the scaling centre of subdomains Ω1 and Ω2 , respectively. Subdomains Ω1 and Ω2 are represented and discretized by the quadratic NURBS ba-

In order to verify the effectiveness of the mortar-based SBIGA in solving electrostatic problems in open domains, as shown in Fig. 15, an electrostatic example in open domain with analytical solution is presented. In the analysis, only half of the open domain is decomposed into bounded domain Ω1 (A–B–C–A) and unbounded domain Ω2 (A–B). O1 and O2 are the scaling centre in subdomains Ω1 and Ω2 , respectively. 27

B. Xue et al.

Engineering Analysis with Boundary Elements 85 (2017) 20–29

Fig. 14. Distribution of potential of U-shaped metal trough.

As shown in Fig. 16, the mortar-based SBIGA provides almost identical distribution of potential compared with that of the analytical solution. Thus, the mortar-based SBIGA can successfully treat the electrostatic problems in open domains without truncating the boundaries or other techniques.

6. Conclusion Fig. 15. Electrostatic problem in open domain.

In this paper, scaled boundary isogeometric analysis is firstly employed to solve the electrostatic problems. Moreover, a mortar-based SBIGA is proposed for treating the electrostatic problems with inhomogeneous media or complex boundaries which cannot satisfy the scaling requirement. The effectiveness, accuracy and convergence properties of SBIGA and the proposed mortar-based SBIGA are verified with benchmark problems in electrostatic field. Compared with SBFEM, SBIGA can exactly describe the boundary, has a rapid convergence rate and provide more accuracy results in solving electrostatic problem with curved

Bounded domain Ω1 and unbounded domain Ω2 are represented and discretized by the quadratic NURBS basis functions independently. The analytical solution of this problem is expressed as:

𝜙=

1 arctan 𝜋

( 𝑥2

2𝑦 + 𝑦2 − 1

) 0 ≤ arctan 𝜑 ≤ 𝜋

(56)

Fig. 16. Distribution of potential in local domain.

28

B. Xue et al.

Engineering Analysis with Boundary Elements 85 (2017) 20–29

boundary. The mortar-based SBIGA, which retain the flexibility of the interface meshes, can successfully solve the electrostatic problems with inhomogeneous media or complex boundaries which cannot satisfy the scaling requirement. A remarkable feature of the mortar-based SBIGA is that, this method is efficient and accurate in solving electrostatic problem in open domains without truncating the boundaries or other techniques. Overall the features of SBIGA and mortar-based SBIGA are shown to make them suitable to analyse other electromagnetic field problems.

[9] Liu J, Lin G. A scaled boundary finite element method applied to electrostatic problems. Eng Anal Boundary Elem 2012;36:1721–32. [10] Wolf JP, Song C. The scaled boundary finite-element method-a primer: derivations. Comput Struct 2000;60:191–210. [11] Song C, Wolf JP. The scaled boundary finite-element method-a primer: solution procedures. Comput Struct 2000;78:211–25. [12] Hughes TJR, Cottrell JA, Bazilevs Y. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Comput Methods Appl Mech Eng 2005;194:4135–95. [13] Anders D, Weinberg K, Reichardt R. Isogeometric analysis of thermal diffusion in binary blends. Comput Mater Sci. 2012;52:182–8. [14] De Lorenzis L, Wriggers P, Hughes TJR. Isogeometric contact: a review. GAMM-Mitteilungen 2014;37:85–123. [15] Espath LFR, Linn RV, Awruch AM. Shape optimization of shell structures based on NURBS description using automatic differentiation. Int J Numer Methods Eng 2011;88:613–36. [16] Buffa A, Sangalli G, Vázquez R. Isogeometric analysis in electromagnetics: B-splines approximation. Comput Methods Appl Mech Eng 2010;199:1143–52. [17] Bernardi C, Maday Y, Patera AT. A new nonconforming approach to domain decomposition: the mortar element method, in Nonlinear Partial Differential Equations and Their Applications. College de France Seminar, Volume XI. Lectures Presented at the Weekly Seminar on Applied Mathematics 1994:13–51. [18] Brivadis E, Buffa A, Wohlmuth B, Wunderlich L. Isogeometric mortar methods. Comput Methods Appl Mech Eng 2015;284:292–319. [19] Seitz A, Farah P, Kremheller J, Wohlmuth BI, Wall WA, Popp A. Isogeometric dual mortar methods for computational contact mechanics. Comput Methods Appl Mech Eng 2016;301:259–80. [20] Wohlmuth BI. Discretization methods and iterative solvers based on domain decomposition. Berlin, Heidelberg, New York: Springer-Verlag; 2001. [21] Lin G, Zhang Y, Hu Z, Zhong H. Scaled boundary isogeometric analysis for 2D elastostatics. Sci China Phys Mech Astron 2014;57:286–300. [22] Natarajan S, Wang J, Song C, Birk C. Isogeometric analysis enhanced by the scaled boundary finite element method. Comput Methods Appl Mech Eng 2015;283:733–62. [23] Piegl L, Tiller W. The NURBS book. 2nd ed. Springer-Verlag New York, Inc.; 1997. [24] Cottrel JA, Hughes TJR, Bazilevs Y. Isogeometric Analysis: towards integration of CAD and CAE. Chichester: Wiley; 2009. [25] Belgacem FB. The mortar finite element method with Lagrange multipliers. Numer Math 1999;84:173–97. [26] Kim J, Youn S. Isogeometric contact analysis using mortar method. Int J Numer Methods Eng 2012;89:1559–81. [27] Fischer KA, Wriggers P, Fischer KA. Mortar based frictional contact formulation for higher order interpolations using the moving friction cone. Comput Methods Appl Mech Eng 2006;195:5020–36.

Acknowledgements This work was supported by the State Key Program of National Natural Science Foundation of China (grant no. 51421064). All supports are gratefully acknowledged. References [1] Abdel-Salam M, El-Mohandes MT. Combined method based on finite differences and charge simulation for calculating electric fields. IEEE Trans Ind Appl 1989;25:1060–6. [2] Dong X, An T. A new FEM approach for open boundary Laplace’s problem. IEEE Trans Microwave Theory Tech 1996;44:157–60. [3] Sumant PS, Cangellaris AC, Aluru NR. A conformal mapping-based approach for fast two-dimensional FEM electrostatic analysis of MEMS devices. Int J Numer Modell Electron Networks Devices Fields 2011;24:194–206. [4] Bo WU. The study and discussion of the electric multipole moment method. J Shangrao Normal Coll 2008;28:20–5. [5] González I, García E, Adana FSD, Cátedra MF. MONURBS: a parallelized fast multipole multilevel code for analyzing complex bodies modeled by NURBS surfaces. Appl Comput Electromag Soc J 2008;23:134–42. [6] Bamji SS, Bulinski AT, Prasad KM. Electric field calculations with the boundary element method. IEEE Trans Electr Insul 1993;28:420–4. [7] He W, Liu Z, Gordon RK, Hutchcraft WE, Yang F, Chang A. A comparison of the element free Galerkin method and the meshless local Petrov-Galerkin method for solving electromagnetic problems. Appl Comput Electromag Soc J 2012;27:620–9. [8] Razmjoo H, Movahhedi M, Hakimi A. An improved truly meshless method based on a new shape function and nodal integration. Int J Numer Modell Electron Networks Devices Fields 2012;25:441–53.

29