Hypersingular integrals

Table of Contents

• Introduction
• Guiggiani's method
  • The original formulation
  • Rong's improvement
• Implementation

Introduction

This tutorial explains how strongly singular and hypersingular integrals are handled in NiHu.

The collocational integral

\( \displaystyle I = \int_{S} K({\bf x}_0, {\bf x}) N({\bf x}) dS_x, \)

is singular if the singular point \( {\bf x}_0 \) is located inside the element domain \( S \).

• If the kernel contains an \( O(1/r^2) \) type singularity in 3D or an \( O(1/r) \) type singularity in 2D, then the integral is strongly singular.
• If the kernel contains an \( O(1/r^3) \) type singularity in 3D or an \( O(1/r^2) \) type singularity in 2D, then the integral is hypersingular.

Guiggiani's method

Guiggiani presented a unified method in 1992 for the accurate numerical evaluation of collocational strongly and hypersingular integrals. Recently, Rong et al published an improved version of the original generic algorithm. NiHu implements the improved version.

The original formulation

As usual, the integration is performed in intrinsic coordinates:

\( \displaystyle I = \int_{\Sigma} K({\bf \xi}_0, {\bf \xi}) N({\bf \xi}) J({\bf \xi}) d\Sigma, \)

where \( \Sigma \) denotes the reference domain of \( S \), \( \xi \) denotes the location vector in the local coordinate system, and \( \xi_0 \) denotes the image of the singular point. \( J(\xi) \) is the Jacobian of the coordinate transform. A further polar coordinate transform is introduced in the local coordinate system around the singular point with the definition

\( \displaystyle \xi = \xi_0 + \rho (\cos \theta, \sin \theta) \)

leading to the double integral

\( \displaystyle I = \int_{0}^{2\pi} \int_{0}^{\bar{\rho}(\theta)} K(\rho, \theta) N(\rho, \theta) J(\rho, \theta) \rho d\rho d \theta = \int_{0}^{2\pi} \int_{0}^{\bar{\rho}(\theta)} F(\rho, \theta) d \rho d \theta \)

The integrand \( F(\rho, \theta)\) is approximated with its Laurent series around the origin of the polar coordinate system (the singular point):

\( \displaystyle F(\rho, \theta) = \frac{F_{-2}(\theta)}{\rho^2} + \frac{F_{-1}(\theta)}{\rho} + O(1) \)

and the truncated expansion is subtracted and added to the integrand to yield

\( \displaystyle I = \int_{0}^{2\pi} \int_{0}^{\bar{\rho}(\theta)} F(\rho, \theta) - \frac{F_{-2}(\theta)}{\rho^2} - \frac{F_{-1}(\theta)}{\rho} d\rho d \theta + \int_{0}^{2\pi} -\frac{F_{-2}(\theta)}{\bar{\rho}(\theta)} + F_{-1}(\theta) \ln |\bar{\rho}(\theta)| d \theta \\ \displaystyle I = \int_{0}^{2\pi} \int_{0}^{\bar{\rho}(\theta)} O(1) d\rho d \theta + \int_{0}^{2\pi} -\frac{F_{-2}(\theta)}{\bar{\rho}(\theta)} + F_{-1}(\theta) \ln |\bar{\rho}(\theta)| d \theta \)

In the last expression both the surface and the line integrals are regular, and can be approximated with standard Gaussian quadrature rules.

The method is fully general, it is valid for any type of curved or distorted surface elements and hypersingular kernels. its only limitation is that the element should be smooth around the singular point. For corners and edges, the line integrals are extended with additional terms, but remain regular.

Rong's improvement

Recently, Rong et al proposed an efficiency improvement to Guiggiani's original formulation. They showed that the intrinsic reference domain \( \Sigma \) can be chosen such that the derivative

\( \displaystyle \left| \lim_{\rho \to 0} \frac{\partial {\bf r}(\rho, \theta)}{\partial \rho} \right| \)

remains constant (independent of \( \theta \)). As this derivative obviously plays an important role in the Laurent coefficients \( F_{-2}(\theta), F_{-1}(\theta) \) of the hypersingular kernel, the improvement yields improved accuracy on highly distorted elements. Furthermore, the new reference domain yields analytical closed form expressions for the angular line integrals.

Implementation

Guiggiani's method with Rong's improvement is implemented in class NiHu::guiggiani in file guiggiani_1992.hpp .

Class NiHu::guiggiani is templated on the trial field (element type and shape function type \( N \)) and the kernel type, as well as on the quadrature orders in radial and tangential directions. The general algorithm can be specialised to a specific kernel type by providing the singularity's Laurent expansion in the class template NiHu::polar_laurent_coeffs.

Note that the template NiHu::polar_laurent_coeffs is not templated on a particular kernel but on the kernel's singularity type. The reason of this differentiation is that different kernels tipically share singularity types. For example, the hypersingular Laplace and Helmholtz kernels differ only in their regular parts. Their singular behaviour, and therefore their Laurent expansion around the singular point are the same.

The class NiHu::polar_laurent_coeffs should implement a static function template called eval that evaluates the laurent coefficients in terms of

• the 2nd order series expansion of the location vector \( {\bf r} \approx {\bf r}_1 \rho + {\bf r}_2 \rho^2 \)
• the 1st order series expansion of the element Jacobian vector \( {\bf J} \approx {\bf J}_0 + {\bf J}_1 \rho \)
• the 1st order series expansion of the shape function \( N \approx N_0 + N_1 \rho \)
• the unit normal vector \( {\bf n} \) at the singular point

These members are obtained from the general guiggiani object taken as parameter by function eval.