Efficient Quadrature Settings for Elliptic PDE’s using a Coupled FEM and BEM Solver in COMSOL Multiphysics

(1)

Efficient Quadrature Settings for Elliptic PDE’s using a Coupled FEM and BEM Solver in COMSOL Multiphysics

G U S T A V F R I B E R G

Master of Science Thesis Stockholm, Sweden 2014

(2)

(3)

Efficient Quadrature Settings for Elliptic PDE’s using a Coupled FEM and BEM Solver in COMSOL Multiphysics

G U S T A V F R I B E R G

Master’s Thesis in Scientific Computing (30 ECTS credits) Master Programme in Scientific Computing (120 credits) Royal Institute of Technology year 2014 Supervisor at KTH was Anna-Karin Tornberg Supervisor at COMSOL AB was Magnus Olsson Examiner was Michael Hanke

TRITA-MAT-E 2014:47 ISRN-KTH/MAT/E--14/47--SE

Royal Institute of Technology School of Engineering Sciences

KTH SCI SE-100 44 Stockholm, Sweden URL: www.kth.se/sci

(4)

(5)

Abstract

By using singular integral kernels based on the fundamental solution, a partial differential equation (PDE) can be rewritten as a boundary integral defined at the boundary of a domain. This requires a linear differential operator with coefficients that are isotropic and homogeneous in space.

In this report, emphasis is put on PDE’s related to elec- tromagnetics i.e., Laplace’s and Helmholtz equation. Both three- and two-dimensional model problems will be inves- tigated.

Galerkin’s method is implemented in order to discretize the domain, now with one less spatial dimension. Hence, the solution is expanded in a series of shape functions ϕ_i whereafter the equation is multiplied by a test function and integrated over the boundary. The resulting matrix ele- ments are double integrals between one shape function ϕ_i and one test function vj integrating vertex i and j in the generated mesh.

Unlike a strict BEM implementation, this report will cover a coupled BEM and FEM solver using Costabel’s Symmetric Coupling. Hence, the resulting system of equa- tions, represented by the stiffness matrix K, consists of both sparse and dense parts originating from the different methods. FEM is usually defined in a domain where there exist non-linearities and BEM is implemented at its boundary in order to simulate an infinite domain as efficiently as possible.

Furthermore, the integrals in K are transformed using two coordinate transforms: one to the reference element and another to avoid the singularity due to the integral kernel. The latter is modified for each case of integration, namely same elements, same edge, same vertex, close elements and distant elements.

The objective of this report is to investigate how the settings for the numerical integration i.e., the quadrature corresponding to the different cases, affect the accuracy of the final solutions to the given PDE’s. However, an ele- ment in K is an integral of a function S which character- istics depend on several things, namely the order of the shape functions, the integral kernel and the element order of the mesh. In order to facilitate the error estimation, the numerical results will be generated from the model problems where the analytical solution is known. An efficient quadrature is achieved when the error originating from the numerical integration of S is small or neglected in compar- ison to the truncation errors i.e., errors originating from meshing and discretization. The thesis is written in close collaboration with the Swedish software company COM-

(6)

(7)

Referat

Genom att använda singulära integralkärnor baserade på den fundamentala lösningen kan en partiell differentialekva- tion (PDE) omskrivas som en randintegral definerad på randen av en domän. Detta kräver en linjär differentialo- perator med tillhörande isotropiska och homogena koffici- enter. Denna rapport fokuserar på PDE:er relaterade till elektromagnetism vilket innebär att fokus kommer läggas på Laplace’s ekvation och Poissons ekvation. Både två- och tre-dimensionella modellproblem kommer att undersökas.

För att diskretisera geometrierna används Galerkin’s metod, nu med en mindre spatiell dimension. Följdaktligen expanderas lösningen i en serie av shapefunktioner ϕ_ivaref- ter ekvationen multipliceras med en testfunktion och sedan integreras över randen. I tre dimensioner blir de resulteran- de matriselementen dubbelintegraler mellan en shapefunc- tion ϕi och en testfunction vj vilka integrerar vertex i och j i den genererade meshen.

Genom att använda Costabel’s Symmetric Coupling kommer denna rapport, till skillnad från en strikt BEM-implementation, behandla en kopplad FEM-/BEM- lösare. Det slutgiltiga systemet av ekvationer, vilket repre- senteras av stiffnessmatrisen K, består därför både av glesa och fyllda delar vilka härstammar från de olika metoderna.

Vanligtvis är FEM definerat i en domän i vilken det exi- sterar olinjäriteter medan BEM är implementerat på do- mänens rand för att så effektivt som möjligt simulera en infinit domän.

Vidare, integralerna i K transformeras genom två ko- ordinattransformer: en till referenselementet och en annan införd i syfte att undvika singulariteten som uppkommer till följd av integralkärnan. Den senare är modifierad för var- je integrationsfall, nämligen samma element, samma kant, samma vertex, närliggande element och avlägsna element.

Målet med rapporten är att undersöka hur inställning- arna för den numeriska integrationen, eller kvadraturen, för de olika integrationsfallen påverkar noggrannheten av den slutgiltiga lösningen till de givna PDE:erna. Ett element i K är en integral av en funktion S vars karaktär beror på flertalet saker, nämligen shapefunktionsordningen, in- tegralkärnan och meshelementens ordning. För att under- lätta feluppskattningen kommer de numeriska resultaten baseras på problem där den analytiska lösningen redan är känd. En effektiv kvadratur är uppnådd när felet som gru- dar sig i den numeriska integrationen av S är litet eller

(8)

ven i nära samarbete med det svenska mjukvaruföretaget COMSOL Multiphysics^® och samtliga numeriska resultat kommer därför genereras med hjälp av denna mjukvara, mer specifikt version 4.4.

(9)

Introduction

The derivations of the mathematical foundations required for creating an integral equation dates back to the 18th century. However, it was not until the middle of the 20th century, when more emphasis was put on developing effective numerical algorithms, the boundary element method (BEM) became more frequently used.

Normally when treating partial differential equations (PDE’s) the problem is defined in a domain Ω and with boundary conditions (BC’s) at the boundary, Γ, of the domain. The boundary conditions defined at Γ could for instance be Dirichlet, Neumann or mixed. By using the linear differential operator L one can define the following PDE using for instance Dirichlet boundary conditions.

( L(u(¯x)) = g x ∈ Ω¯

u(¯x) = f x ∈ Γ¯ (1.0.1)

In order to be able to reformulate such PDE as a boundary integral equation, first one has to compute the fundamental solution G defined as

L(G(¯x)) = δ(¯x) (1.0.2)

or

L(G(¯x, ¯y)) = δ(¯x − ¯y). (1.0.3) The fundamental solution G will in this report be referred to as the Integral Kernel and is different for different choices of L. By integrating both sides of a given PDE over the domain Ω and applying Gauss’s divergence theorem, according to [2], one can rewrite the solution in ¯y as

u(¯y) = Z

S

u(¯x)dG(¯x, ¯y)

dn − G(¯x, ¯y)du(¯x) dn

dS(¯x) (1.0.4) depending on the integral kernel G. By rewriting the PDE as a boundary integral over Γ the solution can be generated infinitely far away by just discretizing the boundary.

(12)

Equation (1.0.4) has too many unknowns in order to be solved in one step.

Equation (1.0.4) is therefore first considered for ¯y ∈ Γ at which it becomes singular when ¯y → ¯x. However, by taking this limit, equation (1.0.4) can be written as

u(¯y) = Z

S

u(¯x)dG(¯x, ¯y)

dn_x − G(¯x, ¯y)du(¯x) dn_x

dS(¯x) ± u(¯y)

2 , y ∈ Γ¯ (1.0.5) where ± denotes from which side of Γ ¯y is approaching, which will be explained in chapter 2.1.4. For a PDE with Dirichlet or Neumann BC’s, one of the parameters in equation (1.0.5) is known, namely u or ^du_dn. By inserting one of these, it is possible to solve the equation at the boundary for the other unknown parameter. Given Dirichlet BC’s, (1.0.5) is solved for ^du_dn. Given Neumann BC’s, (1.0.5) is solved for u. This was the first step. Now, by knowing all the parameters, (1.0.4) can be evaluated for all ¯y. However, to be able to perform step one numerically, the boundary needs to be discretized.

The unknown parameter in step one is discretized by using Galerkin’s method with same test functions as shape functions, meaning this parameter is expanded with respect to a finite number of shape functions taking the value one in each degree of freedom (DOF) or discretized point. The expansion is inserted in equation (1.0.5) and solved for the missing parameter. Furthermore the equation is multiplied by a test function and integrated over the domain Γ resulting in a double integral:

one from the integral formulation and another originating from the multiplication and integration with the test function. Multiplication with each test function yields one additional equation in the system of equations, hence the number of equations are equal to the number of test functions. The resulting matrix K, which often is referred to as the stiffness matrix, is, unlike to FEM, dense where every component is a double integral. In this report, the objective is to find an efficient setting to numerically evaluate these integrals.

This might seem simple but the integral kernels are singular as ¯y → ¯x, hence the double integrals in the stiffness matrix will become difficult to evaluate numerically.

This report will treat the numerical integration of the components in this matrix and how to set up an efficient quadrature depending on for instance shape function order and the integral kernel.

Due to the required linear differential operator L, not all problems can be for- mulated as a boundary integral. However, two solution methods could be coupled. For domains assigned non-linear PDE’s one can implement an expensive FEM based solution. Furthermore, at the boundary of the same domain one could couple FEM with BEM, efficiently simulating an infinite domain. For this COM- SOL Multiphysics^® have implemented an algorithm called Costabel’s Symmetric Coupling [3]. Unsimilar to what was said earlier concerning only one unknown parameter at the boundary, in order to obtain a higher order of generalizability for different BC’s, both parameters u and ^du_dn are in the coupled model considered unknown. One extra unknown requires one more equation, therefore the limit of the flux ^du_dn is introduced according to [3]

(13)

du(¯y)

dn = φ(¯y) = − Z

S

"

dG(¯x, ¯y) dnx

φ(¯x) + d dny

dG(¯x, ¯y) dnx

u(¯x)

#

dS(¯x)±φ(¯y)

2 , y ∈ Γ.¯ (1.0.6) FEM is formulated and discretized inside the domain and at its boundary. BEM is formulated only at the boundary, however not necessarily with the same discretization. With different discretizations follow different test functions, one set of test functions for FEM and another set of test functions for BEM. Coupling is basically achieved by inserting equation (1.0.6) in the formulation for FEM. Hence, the so- lution vector in the resulting coupled system of equations contains both u and ^du_dn and the stiffness matrix K will both have sparse and dense parts originating from the different solution methods, namely FEM and BEM respectively.

(14)

(15)

Chapter 2

BEM and FEM Method

2.1 BEM Implementation

2.1.1 Fundamental Solution to Laplace’s Equation Consider Laplace’s equation

−∆u(¯x) = 0 x ∈ Ω.¯

u(¯x) = f x ∈ Γ.¯ (2.1.1)

In order to successfully rewrite a PDE as an Integral equation one first has to find its fundamental solution G defined as

L(u) = δ. (2.1.2)

For Laplace’s equation

−∆_xG(¯x, ¯y) = δ(¯x − ¯y) (2.1.3) where δ(¯x − ¯y) denotes the Dirac delta function and represents a unit source term at the point ¯y. Integrating both sides of equation (2.1.3) over the domain V yields

Z

V

−∆G(¯x, ¯y) dV = Z

V

δ(¯x − ¯y) dV = 1. (2.1.4) Gauss’s divergence theorem results in the following where dS is the surface element of a artificial sphere surrounding the source term. Due to the symmetric source term, the solution is hence symmetric resulting in r dependence where r is defined as the euclidean norm of the integrating vector.

(16)

1 = Z

V

−∆G(¯x, ¯y) dV = Z

V

−∇ · ∇G(¯x, ¯y) dV = Z

S

−∇G(¯x, ¯y) · ˆn dS =

r dependence

(2.1.5)

= Z

S

−dG dr dS =

Sphere of radius a

= −4πa²dG dr r=a

(2.1.6) which implies

dG

dr = − 1

4πr² ⇒ G(r) = 1

4πr + C = 1

4π|¯x − ¯y| (2.1.7) on a sphere with radius r centered around the source point. The constant C is by physical reasons put to zero since the potential U should approach zero as ¯x → ∞, which defines the boundary condition. As one can deduce from equation (2.1.7) the fundamental solution, also referred to as the Integral Kernel, is singular as ¯x approaches ¯y. This whole report will revolve around this property. Equation (2.1.7) is valid for ¯x ∈ R³. For ¯x ∈ R² the fundamental solution takes the form

u(r) = −log(r)

2π . (2.1.8)

2.1.2 Integral Equation of Laplace’s Equation

Applying

∆ab = ∇(∇ab) − ∇a∇b (2.1.9)

and by combining the fundamental solution (2.1.3) with the solution to (2.1.1) one obatin

Z

V

u(¯x)(−∆xG(¯x, ¯y))dV = Z

V

[−∇_x(u∇_xG) + ∇xu∇xG] dV. (2.1.10) Z

V

G(¯x, ¯y)(−∆xu(¯x))dV = Z

V

[−∇_x(G∇_xu) + ∇xG∇xu] dV. (2.1.11) Subtracting (2.1.11) from (2.1.10) gives us:

Z

V





−u(¯x) ∆G(¯x, ¯y)

| {z }

= −δ(¯x − ¯y)

+G(¯x, ¯y) ∆u(¯x)

| {z }

=0





dV = Z

V

[−∇(u∇G) + ∇(G∇u)] dV =

=

Gauss’s divergence theorem

=

= Z

S

G(¯x, ¯y)du(¯x)

dn − u(¯x)dG(¯x, ¯y) dn

dS(¯x).

(17)

2.1. BEM IMPLEMENTATION

According to (2.1.1) and (2.1.3), ∆u = 0 and ∆G = −δ(¯x − ¯y), hence the resulting equation becomes

u(¯y) = Z

S

G(¯x, ¯y)du(¯x)

dS(¯x). (2.1.12) According to [2], Table (2.1.13) shows the two- and three-dimensional integral ker- nels G for Laplace’s equation.

Spatial Dimension Integral Kernel G 2 −^log(|¯_2π^x−¯^y|) 3 _4π|¯_x−¯¹ _y|

. (2.1.13)

2.1.3 Integral Equation of Helmholtz Equation

For Helmholtz equation, just switch the differential operator L defined as

L(u) = (−∆ − k²)u. (2.1.14)

The boundary condition for Helmholtz equation are also set infinitely far away from the source term and is based on the energy flux defined by the Pyonting vector. The energy flux should be conserved and proportional to _r¹₂. In three dimensions, this implies that the solution u is proportional to ¹_r. Furthermore, in two dimensions, u should be proportional to ^√¹

r as r → ∞.

Z

V

u(¯x)

− ∆G(¯x, ¯y) − k²G(¯x, ¯y)

dV = Z

V

h−∇(u∇G) + ∇u∇G − uk²GⁱdV.

Z

V

G(¯x, ¯y)

− ∆u(¯x) − k²u(¯x, ¯y)

dV = Z

V

h−∇(G∇u) + ∇G∇u − Gk²uⁱdV.

(2.1.15) Subtraction results in cancellation of ∇u∇G and uk²G.

Z

V





 u(¯x)

− ∆G(¯x, ¯y) − k²G(¯x, ¯y)

| {z }

=δ(¯x,¯y)

−G(¯x, ¯y)

− ∆u(¯x) − k²u(¯x, ¯y)

| {z }

=0





 dV

= Z

V

[−∇(u∇G) + ∇(G∇u)] dV =

Gauss’s divergence theorem

=

= Z

S

G(¯x, ¯y)du(¯x)

dS(¯x).

(18)

In analogy with the derivation from Laplace’s equation we obtain the following equation defined at the boundary.

u(¯y) = Z

S

G(¯x, ¯y)du(¯x)

dS(¯x). (2.1.16) According to equations (2.1.12) and (2.1.16), what differs between the integral formulations for Laplace’s equation and Helmholtz equation is just the integral kernel G. According to [2], the two- and three-dimensional integral kernels for Helmholtz equation are shown in Table (2.1.17).

Spatial Dimension Integral Kernel G 2 ₄ⁱH₀⁽¹⁾(k|¯x − ¯y|) 3 ^e_4π|¯^ik|¯_x−¯^{x− ¯}^y|_y|

. (2.1.17)

Here H₀⁽¹⁾(x) denotes a Hankel function from which also the analytical solution is derived through convolution with the source function.

2.1.4 BEM Formulation at the Boundary

According to section 1, in order to have a general solver that easily can manage different boundary conditions, both u and φ := ^du_dn are always considered unknown at the boundary. However, this implies two unknowns at the boundary which requires one more equation. The normal derivative of (2.1.12) is therefore taken resulting in,

φ = ∂u(¯y)

∂n = d dn

Z

S

u(¯x)dG(¯x, ¯y)

dS(¯x), y ∈ Ω.¯ (2.1.18) The following derivation shows the limit for equation (1.0.4) as ¯x approaches ¯y at the boundary. The integral is split up into two surface integrals: one over the non-singular region and another for the singular region. In the singular region a hemisphere with radius is placed around the singularity giving us below limit for

 → 0. The normal of the hemisphere is |¯x− ¯y| which coincides with the denominator in the kernel function G. The radius of the hemisphere is denoted and is a scalar.

Depending on the orientation of the hemisphere i.e., if ¯x is approaching from the inside or the outside there will be a change in sign due to the limits of integration for θ. Starting with equation (2.1.12) gives us

Z

S

u(¯x)dG(¯x, ¯y)

dS(¯x) =

→0lim Z

S−S

u(¯x)dG(¯x, ¯y)

dS(¯x)+

Z

S

u(¯x)dG(¯x, ¯y)

dS(¯x) (2.1.19)

(19)

2.1. BEM IMPLEMENTATION

With

G(¯x, ¯y) = − 1

4π|¯x − ¯y| (2.1.20)

the second integral on the right hand side of equation (2.1.19) becomes

lim→0

Z π/2 θ=0

Z 2π φ=0

u(¯x)dG(¯x, ¯y)

²sin(θ)dθdφ = lim→0

Z π/2 θ=0

Z 2π φ=0

u(¯x)d

d

− 1

4π|¯x − ¯y|

+ 1

4π|¯x − ¯y|

du(¯x) d

²sin(θ)dθdφ =

|¯x − ¯y| = 

= lim→0

Z π/2 θ=0

Z 2π φ=0

u(¯x)d

d

− 1 4π

+ 1

4π

du(¯x) d

Z π/2 θ=0

Z 2π φ=0

u(¯x) 1

4π² + 1 4π

du(¯x) d

Z π/2 θ=0

Z 2π φ=0

u(¯x) 1

4π + 4π

du(¯x) d

| {z }

→0

sin(θ)dθdφ =

lim→0

Z π/2 θ=0

Z 2π φ=0

u(¯x) 1

4π sin(θ)dθdφ =

u = u0+ (u − u₀)

= lim→0

Z π/2 θ=0

Z 2π φ=0

u₀(¯x) + u(¯x) − u₀(¯x)

4π sin(θ)dθdφ =

lim→0

Z π/2 θ=0

Z 2π φ=0

u0(¯x)

4π sin(θ)dθdφ + lim

→0

Z π/2 θ=0

Z 2π φ=0

u(¯x) − u0(¯x)

4π sin(θ)dθdφ

| {z }

u0(¯x)→u(¯x) as →0

=

lim

→0

Z π/2 θ=0

Z 2π φ=0

u₀(¯x)

4π sin(θ)dθdφ = lim

→0

2πu₀(¯y)

4π = ±u₀(¯y) 2

(2.1.21) This gives us the following result

u(¯y) = Z

S

u(¯x)dG(¯x, ¯y)

dn_x − G(¯x, ¯y)du(¯x) dn_x

dS(¯x) ± u₀(¯y)

2 , y ∈ Γ.¯ (2.1.22) In analogy with the limit for (2.1.12), for y ∈ Γ, (2.1.18) becomes

φ(¯y) = − Z

S

"

dG(¯x, ¯y) dnx

φ(¯x) + d dny

dG(¯x, ¯y) dnx

u(¯x)

#

dS(¯x) ±φ(¯y)

2 , y ∈ Γ¯ (2.1.23) where ± denotes from which side of Γ ¯y is approaching.

(20)

2.2 Coordinate Transform to Reference Triangle

2.2.1 The Reference Triangle

In three dimensions, all integrals are mapped to the reference triangle which has its corners located in (0,0), (1,0) and (1,1) according to Figure 2.1.

(0,0) (1,0)

(1,1)

Figure 2.1: All integrals are mapped to the reference triangle shown in this Figure.

In some books (1,1) is moved to (0,1).

The transform depends on the location of the surface panel but also its shape. Linear elements is the simplest possible case but another type of mapping, often referred to as isoparametric, has to be implemented if the edges of the global element are curved.

2.2.2 Linear Surface Panels

In a linear transform the reference triangle maps to a surface panel in the global coordinate system that has straight edges. Hence, using a linear transform it is not possible to map a geometry with curved edges. The coordinate transform should be chosen carefully in order to reduce numerical errors to minimal cost. Using linear basis functions

x = ¯¯ a(1 − ˆx₁) + ¯b(ˆx₂) + ¯c(ˆx₁− ˆx₂) (2.2.1)

where ¯a, ¯b and ¯c are the corners of a surface panel in the global system.

(21)

2.2. COORDINATE TRANSFORM TO REFERENCE TRIANGLE

(0,0) (1,0)

(1,1)

a b

c

Figure 2.2: Linear mapping of the reference triangle.

The Jacobian J is given by

J = d¯x dˆx =







−a_x+ c_x bx− c_x

−a_y+ c_y b_y− c_y

−a_z+ c_z b_z− c_z





, (2.2.2)

where the determinant is defined as

det(J ) =^q(J_1,1J_2,2− J_2,1J_1,2)²+ (J_2,1J_3,2− J_3,1J_2,2)²+ (J_3,1J_1,2− J_1,1J_3,2)². (2.2.3) In the linear case (2.2.2) becomes constant and does not affect the method of integra- tion for the elements in K. However, when using curved surface panels the mapping cannot be linear, hence the determinant of the Jacobian becomes dependent on the space variables ˆx, ˆy.

2.2.3 Quadratic Surface Panels

For surface panels with curved edges one has to interpret a non-linear transform, which is an example of an isoparametric mapping. For a surface panel with three vertices and curved edges with one interpolating point on each edge we obtain the following coordinate transform. Consider it as six basis functions that all equal 1 at one point in the triangle and 0 at all other points.

(0,0) (1/2,0)

(1/2,1/2) (1,1/2) (1,1)

(1,0)

a b

c

d e

f

Figure 2.3: Quadratic mapping of reference triangle.

(22)

x =2¯¯ a(1 − ˆx₁)(1

2 − ˆx₁) + 2¯bˆx₂(ˆx₂−1

2) + ¯c(ˆx₁− ˆx₂)(ˆx₁− ˆx₂−1 2)+

4 ¯dˆx₂(1 − ˆx₁) + 2¯eˆx₂(ˆx₁− ˆx₂) + 4 ¯f (ˆx₁− ˆx₂)(1 − ˆx₁)

(2.2.4)

The resulting elements of the Jacobian are listed below.

J_1,1= 2¯a_x−3

2 + 2ˆx₁+ ¯c_x2ˆx₁− 2ˆx₂− 1 2

+ 4 ¯d_x− ˆx₂+ 2¯e_xxˆ₂+ 4 ¯f_x1 − 2ˆx₁+ ˆx₂ (2.2.5) J_1,2= 2¯b_x2ˆx₂−1

2

+ ¯c_x− 2ˆx₁− 2ˆx₂−1 2

+ 4 ¯d_x1 − ˆx₁+ 2¯e_xxˆ₁− 2ˆx₂+ 4 ¯f_x− 1 − ˆx₁ (2.2.6)

J2,1= 2¯ay

−3

2 + 2ˆx1 + ¯cy

2ˆx1− 2ˆx2−1 2

+ 4 ¯dy

− ˆx2

+ 2¯eyxˆ2+ 4 ¯fy

1 − 2ˆx1+ ˆx2 (2.2.7) J2,2= 2¯by

2ˆx2−1 2

+ ¯cy

− 2ˆx1− 2ˆx2−1 2

+ 4 ¯dy

1 − ˆx1 + 2¯ey

xˆ1− 2ˆx2 + 4 ¯fy

− 1 − ˆx1 (2.2.8)

J_3,1= 2¯a_z−3

2+ 2ˆx₁+ ¯c_z2ˆx₁− 2ˆx₂−1 2

+ 4 ¯d_z− ˆx₂+ 2¯e_zxˆ₂+ 4 ¯f_z1 − 2ˆx₁+ ˆx₂ (2.2.9) J3,2= 2¯bz

2ˆx2−1 2

+ ¯cz

− 2ˆx1− 2ˆx2−1 2

+ 4 ¯dz 1 − ˆx1

+ 2¯ez

xˆ1− 2ˆx2 + 4 ¯fz

− 1 − ˆx1 (2.2.10)

For the non-linear transform the determinant is defined in accordance with (2.2.3).

To create an even more accurate transform i.e., mesh elements of higher order than two, one simply adds more points to the transform. The method is often referred to as Lagrange interpolation on a triangle. Table (2.2.11) shows the number of interpolating points and shape functions required for one mesh element in order to obtain a certain geometric order.

Linear Quadratic Cubic . . .

Interpolating points 3 6 9 . . .

Number of Basis functions 3 6 9 . . .

(2.2.11)

2.3 Shape Functions

With the defined basis for the transform one has to derive the shape functions ϕ representing the solution. Depending on requested accuracy and problem character,

(23)

2.3. SHAPE FUNCTIONS

different shape functions are used. Common for all types of shape functions is the following restriction which holds for Lagrangian shape function elements

ϕ_i(¯x_j) = δ_ij (2.3.1)

showing that every shape function only takes the value 1 at just one vertex in the mesh.

2.3.1 Linear Shape Functions

The following example shows how to determine linear shape functions for triangular surface panels. Begin with a linear ansatz

σ(ˆx₁, ˆx₂) = c₁+ c₂xˆ₁+ c₃xˆ₂. (2.3.2)

One corner in the reference triangle maps to one corner in the global triangle.

σ(0, 0) = c1 = σ₁ (2.3.3)

σ(1, 1) = c₁+ c₂+ c₃= σ₂ (2.3.4) σ(1, 0) = c₁+ c₂= σ₃ (2.3.5)

Three equations, three unknowns and reordering of the terms gives us

σ(ˆx1, ˆx2) = σ₁(1 − ˆx1)

| {z }

ϕ1(ˆx)

+σ₂(ˆx2)

| {z }

ϕ2(ˆx)

+σ₃(ˆx1− ˆx2)

| {z }

ϕ3(ˆx)

. (2.3.6)

Note the similarities with the basis functions for the coordinate transform. As one can deduce from (2.3.6) equation (2.3.1) is valid meaning

(0,0) (1,1) (1,0)

ϕ1 1 0 0

ϕ2 0 1 0

ϕ₃ 0 0 1

(24)

Figure 2.4: Linear Shape Functions plotted on the Reference Triangle.

2.3.2 Quadratic Shape Functions Begin with a quadratic ansatz.

σ(ˆx1, ˆx2) = c₁+ c₂xˆ1+ c₃xˆ2+ c₄xˆ²₁+ c₅xˆ²₂+ c₆xˆ1xˆ2. (2.3.7)

Six unknowns require six equations. Since the triangle only has 3 corners, we introduce 3 new points located at the middle of each edge. Still the resulting shape functions should fulfill (2.3.1).

σ(0, 0) = c₁= σ₁ (2.3.8)

σ(1, 1) = c₁+ c₂+ c₃+ c₄+ c₅+ c₆= σ₂ (2.3.9) σ(1, 0) = c1+ c₂+ c₄ = σ₃ (2.3.10) σ(1

2,1

2) = c₁+c₂ 2 +c₃

2 +c₄ 4 +c₅

4 +c₆

4 = σ₄ (2.3.11)

σ(1,1

2) = c₁+ c₂+c3

2 + c₄+c5

4 + c6

2 = σ₅ (2.3.12)

σ(1

2, 0) = c1+c₂ 2 +c₄

2 = σ₆ (2.3.13)

The resulting function on that triangle follows

(25)

2.4. GALERKIN’S METHOD

σ(ˆx1, ˆx2) =σ₁2(1 − ˆx1)(1 2 − ˆx1)

| {z }

ϕ1(ˆx)

+σ₂2ˆx2(ˆx2−1 2)

| {z }

ϕ2(ˆx)

+

σ3(ˆx1− ˆx2)(ˆx1− ˆx2−1 2)

| {z }

ϕ3(ˆx)

+σ₄4ˆx2(1 − ˆx1)

| {z }

ϕ4(ˆx)

+

σ₅2ˆx₁(ˆx₁− ˆx₂)

| {z }

ϕ5(ˆx)

+σ₆4(ˆx₁− ˆx₂)(1 − ˆx₁)

| {z }

ϕ6(ˆx)

(2.3.14)

According to equation (2.3.1), the following holds.

(0,0) (1,1) (1,0) (¹₂,¹₂) (1, ¹₂) (¹₂, 0)

ϕ₁ 1 0 0 0 0 0

ϕ2 0 1 0 0 0 0

ϕ3 0 0 1 0 0 0

ϕ₄ 0 0 0 1 0 0

ϕ5 0 0 0 0 1 0

ϕ6 0 0 0 0 0 1

Figure 2.5: Quadratic Shape Functions on the Reference Triangle.

Figure 2.6: Quadratic Shape Functions on the Reference Triangle.

2.4 Galerkin’s Method

Equation (2.1.22) and (1.0.6) are solved at the boundary in order to determine the unknown parameters, u and _dn^du on Γ, whereafter equation (2.1.12) gives the solution in a point ¯y. However, the system needs to be discretized. Comsol Multiphysics^® uses a method based on Galerkin’s orthogonality, hence u and φ are expanded as sets of shape functions. Normally the summations of u and φ do not consist of

(26)

the same number of terms. u is defined inside the domain and at the boundary.

With A degreees of freedom (DOF) inside the domain and B DOF at the boundary, the expansion of u is a summation from i = 1 to i = A + B. Furthermore, φ is only defined at the boundary. Discretizing φ with C DOF results in an expansion consisting of C terms. Another discretization also implies different shape functions θ. Equations (2.4.1) and (2.4.2) show the expansions.

u(¯y) =

A+B

X

i=1

uiϕi(¯y). (2.4.1)

du

dn = φ(¯y) =

C

X

k=1

φ_kθ_k(¯y). (2.4.2)

In COMSOL Multiphysics^®, for both formulations i.e., FEM and BEM, the set of test functions are set equal to the set of shape functions.

2.4.1 BEM

Equations (2.4.1) and (2.4.2) are inserted in equation (2.1.22) whereafter both sides of the equation are multiplied by a testfunction ψ_j and integrated over the boundary Γ_y. In three dimensions the result is a double surface integral: one over surface element S(¯x) and another over surface element S(¯y). In two dimensions this is a double integral over two line segments. For three-dimensional problems:

A+B

X

i=1

Z

S(¯y)

ϕ_i(¯x)u_i

2 ψ_j(¯y)dS(¯y) =

A+B

X

i=1

Z

S(¯y)

Z

S(¯x)

dG(¯x, ¯y) dnx

u_iϕ_i(¯x)ψ_j(¯y)dS(¯x)dS(¯y)−

C

X

k=1

Z

S(¯y)

Z

S(¯x)

G(¯x, ¯y)φkθk(¯x)ψj(¯y)dS(¯x)dS(¯y), x, ¯¯ y ∈ Γ.

(2.4.3)

The result is a system of C equations, one for each test function. Unsimilar to the FEM where the matrix becomes sparse, G(¯x, ¯y) connects the two integrals on the right hand side resulting in a full system. Thus, element j, i of such matrix is therefore the integral of vertex i and j in the genrated mesh. However, as j approaches i, these surface integrals become singular. This problem only arises in BEM and emphasis is put on this later in section 2.6. Using following notations for R, B and W

(27)

2.4. GALERKIN’S METHOD

R(u) = Z

S

dG(¯x, ¯y)

dn_x u_iϕ_i(¯x)dS(¯x), (2.4.4) B(φ) =

Z

S

G(¯x, ¯y)φ_kϕ_k(¯x) dS(¯x), (2.4.5) W (u) = − d

dn_y Z

S

dG(¯x, ¯y)

dn_x u_iϕ_i(¯x)dS(¯x), (2.4.6) equation j in (2.4.3) takes the form

u 2, ψj

Γ

−

R(u), ψj

Γ

+

B(φ), ψj

Γ

= 0 (2.4.7)

where u is a summation from 1 → A + B and φ is a summation from 1 → C.

Furthermore, (·, ·) denotes the L₂ scalar product over the specified domain.

2.4.2 FEM Formulation to Laplace’s Equation Given Poisson’s equation

−∆u(¯x) = g(¯x) x ∈ Ω.¯ (2.4.8)

Insert (2.4.1), multiply both sides by another test function v_j(¯y) and integrate over the domain Ω

−

A+B

X

i=1

Z

Ω

u_i∆ϕ_i(¯x)v_j(¯x)dΩ = Z

Ω

g(¯x)v_j(¯x)dΩ. (2.4.9) Applying equation (2.1.9) on equation (2.4.9) results in

−

A+B

X

i=1

Z

Ω

ui∇ ∇ϕ_i(¯x)vj(¯x)− u_i∇ϕ_i(¯x)∇vj(¯x)dΩ = Z

Ω

g(¯x)vj(¯x)dΩ. (2.4.10) Applying Gauss’s divergence theorem results in

−

A+B

X

i=1

Z

S

u_i∇ϕ_i(¯x)v_j(¯x)ˆndS(¯x) +

A+B

X

i=1

Z

Ω

u_i∇ϕ_i(¯x)∇v_j(¯x)dΩ = Z

Ω

g(¯x)v_j(¯x)dΩ.

(2.4.11)

−

C

X

k=1

Z

S

φ_kθ_k(¯x)v_j(¯x)dS(¯x) +

A+B

X

i=1

Z

Ω

u_i∇ϕ_i(¯x)∇v_j(¯x)dΩ = Z

Ω

g(¯x)v_j(¯x)dΩ.

(2.4.12) With L₂ norm notation equation j in the system can be written as

∇u, ∇v_j

Ω

−

φ, vj

Γ

=

g, vj

Ω

. (2.4.13)

(28)

2.4.3 FEM Formulation to Helmholtz Equation

−∆u(¯x) − k²u = g(¯x) x ∈ Ω¯ (2.4.14) In accordance with section 2.4.2, insert equation (2.4.1), multiply both sides by a test function v_j(¯y) and integrate over the domain Ω

−

A+B

X

i=1

Z

Ω

hu_i∆ϕ_i(¯x)v_j(¯x) + k²u_iϕ_i(¯x)v_j(¯x)ⁱdΩ = Z

Ω

g(¯x)v_j(¯x)dΩ. (2.4.15)

Applying (2.1.9) gives

−

A+B

X

i=1

Z

Ω

h

ui∇ ∇ϕ_i(¯x)vj(¯x)− u_i∇ϕ_i(¯x)∇vj(¯x)v + k²uiϕi(¯x)vj(¯x)ⁱdΩ =

= Z

Ω

g(¯x)v_j(¯x)dΩ.

(2.4.16) Gauss’s divergence theorem.

−

A+B

X

i=1

Z

S

ui∇ϕ_i(¯x)vj(¯x)ˆndS(¯x) +

A+B

X

i=1

Z

Ω

ui∇ϕ_i(¯x)∇vj(¯x)−

k²uiϕi(¯x)vj(¯x)dΩ = Z

Ω

g(¯x)vj(¯x)dΩ.

(2.4.17)

−

C

X

k=1

Z

S

φ_kθ_k(¯x)v_j(¯x)dS(¯x) +

A+B

X

i=1

Z

Ω

hu_i∇ϕ_i(¯x)∇v_j(¯x) − k²u_iϕ_i(¯x)v_j(¯x)ⁱdΩ =

= Z

Ω

g(¯x)v_j(¯x)dΩ.

(2.4.18) With L₂ norm notation equation j in the system can be written as

∇u, ∇v_j

Ω

−

φ, v_j

Γ

− k²

u, v_j

Ω

=

g, v_j

Ω

. (2.4.19)

2.5 Costabel’s Symmetric Coupling

Consider two domains Ω⁻ and Ω⁺, where (−) denotes the interior domain and (+) denotes the exterior domain. When treating electromagnetic problems, the field or the potential far away from the source is often of interest. Since FEM requires a fully discretized domain this method therefore becomes extremely expensive for these

(29)

2.5. COSTABEL’S SYMMETRIC COUPLING

kinds of problems, hence a model that couples FEM and BEM becomes relevant.

FEM solves the inhomogeneous part of the PDE in the interior domain Ω⁻ and BEM is applied at the boundary in order to reach infinity as effectively as possible in the exterior domain Ω⁺. Figure 2.7 shows a typical problem where a combination of these two methods is applicable.

u = g, x ∊ ⁻ u = 0, x ∊ ⁺ u = f, x ∊ L

L

Figure 2.7

According to Figure 2.7, one preferably chooses to use FEM in Ω⁻and BEM in Ω⁺. However, these two formulations have to be solved in one system of equations since the solutions depend on each other. Hence, equation (2.4.7) and (2.4.13) have to be coupled. COMSOL Multiphysics^® uses a method called Costabel’s Symmetric Coupling. After discretizing the equations using Galerkin’s method one can write the PDE as a matrix-vector multiplication where the matrix often is referred to as the stiffness matrix. In this report the stiffness matrix is denoted K. Source terms, g, originating from FEM, located in the interior domain Ω⁻, are defined at the right hand side in the system of equations. The following matrix assembly of K is defined for Laplace’s equation. For Helmholtz equation one simply switch the integral kernel according to (2.1.17) and add the extra term k²(u, v_j)_Ω from equation (2.4.19). This extra non-singular term is related to FEM and is therefore located in the sparsity pattern of K for the first A + B equations. According to [3], rewriting (1.0.6) in terms of the R and W gives us

φ = −R(φ) − W (u) +φ

2, y ∈ Γ.¯ (2.5.1)

According to equation (2.1.21), the plus sign in (2.5.1) is due to the integration of θ. The solution is sought in Ω⁺, hence the limit is defined when ¯y → ¯x from the outside. Inserting (2.5.1) in (2.4.13)

∇u, ∇v_j

Ω

+

R(φ) − φ 2, vj

Γ

+

W (u), v

Γ

=

g, vj

Ω

(2.5.2)

where again u and φ are summations. The following two equations is the result of the Costabel’s coupling. Below states (2.4.3) one more time on integral form.

Efficient Quadrature Settings for Elliptic PDE’s using a Coupled FEM and BEM Solver in COMSOL Multiphysics

Efficient Quadrature Settings for Elliptic PDE’s using a Coupled FEM and BEM Solver in COMSOL Multiphysics

Efficient Quadrature Settings for Elliptic PDE’s using a Coupled FEM and BEM Solver in COMSOL Multiphysics

Abstract

Referat

Contents

Chapter 1

Introduction

Chapter 2

BEM and FEM Method

2.1 BEM Implementation

2.2 Coordinate Transform to Reference Triangle

2.3 Shape Functions

2.4 Galerkin’s Method

2.5 Costabel’s Symmetric Coupling