Integral equation formulation of the biharmonic Dirichlet problem

TL;DR

This paper introduces a new integral representation for solving the biharmonic Dirichlet problem, improving robustness and accuracy, especially on high-curvature domains, by leveraging a connection to Stokes flow.

Contribution

The authors develop a novel integral formulation for the biharmonic Dirichlet problem using a Stokes problem analogy, applicable to complex domains, enhancing computational stability.

Findings

New integral representation with better kernel behavior on high-curvature domains

Augmentation techniques for handling simply and multiply connected domains

Numerical examples demonstrating improved robustness and accuracy

Abstract

We present a novel integral representation for the biharmonic Dirichlet problem. To obtain the representation, the Dirichlet problem is first converted into a related Stokes problem for which the Sherman-Lauricella integral representation can be used. Not all potentials for the Dirichlet problem correspond to a potential for Stokes flow, and vice-versa, but we show that the integral representation can be augmented and modified to handle either simply or multiply connected domains. The resulting integral representation has a kernel which behaves better on domains with high curvature than existing representations. Thus, this representation results in more robust computational methods for the solution of the Dirichlet problem of the biharmonic equation and we demonstrate this with several numerical examples.

Figures

Equations

$${\Delta }^{2}w$$

$$w$$

$$\frac{\partial w}{\partial n}$$

$${\Delta }^{2}w$$

$$\frac{\partial w}{\partial \tau }$$

$$\frac{\partial w}{\partial n}$$

$$-{\Delta }^{2}w=q$$

$$w=0$$

$$\frac{\partial w}{\partial n}=0$$

$$\Delta u$$

$$u$$

$$v(\mathbf{x})=-\frac{1}{2\pi }{\int }_{\Omega }\log |\mathbf{x}-\mathbf{y}|\tilde{f}(\mathbf{y})dy,$$

$${\Delta }^{2}w=0$$

$$w=f$$

$$\frac{\partial w}{\partial n}=\frac{\partial f}{\partial n}$$

$$w=0$$

$$\frac{\partial w}{\partial n}=0$$

$$w(\mathbf{x})={\int }_{\Gamma }{K}_1^F(\mathbf{x},\mathbf{y}){\sigma }_1(\mathbf{y})dS(\mathbf{y})+{\int }_{\Gamma }{K}_2^F(\mathbf{x},\mathbf{y}){\sigma }_2(\mathbf{y})dS(\mathbf{y}),$$

$$G^B(\mathbf{x},\mathbf{y})=\frac{1}{8\pi }|\mathbf{x}-\mathbf{y}{|}^2\log |\mathbf{x}-\mathbf{y}|.$$

$$K_1^F(\mathbf{x},\mathbf{y})$$

$$K_2^F(\mathbf{x},\mathbf{y})$$

$$K_1^F(\mathbf{x},\mathbf{y})$$

$$K_2^F(\mathbf{x},\mathbf{y})$$

$$\begin{array}{rl}\left(\begin{array}{c}f(\mathbf{x})\\ g(\mathbf{x})\end{array}\right)& ={\int }_{\Gamma }\left(\begin{array}{cc}{K}_{11}^F(\mathbf{x},\mathbf{y})& K_{12}^F(\mathbf{x},\mathbf{y})\\ K_{21}^F(\mathbf{x},\mathbf{y})& K_{22}^F(\mathbf{x},\mathbf{y})\end{array}\right)\left(\begin{array}{c}{\sigma }_1(\mathbf{y})\\ {\sigma }_2(\mathbf{y})\end{array}\right)dS(\mathbf{y})\\ & +\left(\begin{array}{cc}1/2& 0\\ -\kappa (\mathbf{x})& 1/2\end{array}\right)\left(\begin{array}{c}{\sigma }_1(\mathbf{x})\\ {\sigma }_2(\mathbf{x})\end{array}\right),\end{array}$$

$$K_{21}^F(\mathbf{x},\mathbf{y})$$

$$=\frac{1}{\pi }\left(-3\frac{(\mathbf{r}\cdot \mathbf{n}(\mathbf{y}){)}^2(\mathbf{n}(\mathbf{x})\cdot \mathbf{n}(\mathbf{y}))}{(\mathbf{r}\cdot \mathbf{r}{)}^2}+4\frac{(\mathbf{r}\cdot \mathbf{n}(\mathbf{y}){)}^3(\mathbf{r}\cdot \mathbf{n}(\mathbf{x}))}{(\mathbf{r}\cdot \mathbf{r}{)}^3}\right),$$

$$K_{22}^F(\mathbf{x},\mathbf{y})$$

$$=\frac{1}{\pi }\left(\frac{(\mathbf{r}\cdot \mathbf{n}(\mathbf{y}))(\mathbf{r}\cdot \mathbf{n}(\mathbf{x}))}{\mathbf{r}\cdot \mathbf{r}}-\frac{(\mathbf{r}\cdot \mathbf{n}(\mathbf{y}){)}^2(\mathbf{r}\cdot \mathbf{n}(\mathbf{x}))}{(\mathbf{r}\cdot \mathbf{r}{)}^2}\right).$$

$$-\Delta \mathbf{u}+\nabla p$$

$$\nabla \cdot \mathbf{u}$$

$$\mathbf{u}$$

$${\boldsymbol{u}}=\nabla^{\perp}w=\left(\begin{array}[]{r}\frac{\partial w}{\partial x_{2}}\\ -\frac{\partial w}{\partial x_{1}}\end{array}\right),$$

$$w\left(x_1,x_2\right)=\text{Re}\left(\bar{z}\varphi \left(z\right)+\chi \left(z\right)\right),$$

$$\frac{\partial w}{\partial x_1}+i\frac{\partial w}{\partial x_2}=\varphi \left(z\right)+z\overline{{\varphi }^{{}^{\prime }}\left(z\right)}+\overline{\psi \left(z\right)}.$$

$$G_{ij}\left(\mathbf{x},\mathbf{y}\right)=\frac{1}{4\pi }\left[-\log \left|\mathbf{x}-\mathbf{y}\right|{\delta }_{ij}+\frac{\left(x_i-y_i\right)\left(x_j-y_j\right)}{{\left|\mathbf{x}-\mathbf{y}\right|}^2}\right]i,j\in \{1,2\}.$$

$${\left[S_{\Gamma }\mu \right]}_i\left(\mathbf{x}\right)={\int }_{\Gamma }{G}_{ij}\left(\mathbf{x},\mathbf{y}\right){\mu }_j\left(\mathbf{y}\right)dS(\mathbf{y})i=1,2.$$

$${\varphi }_S(z)$$

$${\chi }_S(z)$$

$$+\frac{1}{8\pi }{\int }_{\Gamma }\overline{\xi }\rho (\xi )\log \left(\xi -z\right)dS(\xi ),$$

$${\psi }_S(z)$$

$$\begin{array}{rl}{w}_S(z)& =\text{mbox}Re[\frac{1}{4\pi }{\int }_{\Gamma }\text{mbox}Re\left[\overline{\left(\xi -z\right)}\rho \left(\xi \right)\right]\log \left(\xi -z\right)dS(\xi )\\ & -\frac{1}{8\pi }{\int }_{\Gamma }\overline{\rho \left(\xi \right)}\left(\xi -z\right)dS(\xi )+\overline{z}\frac{1}{8\pi }{\int }_{\Gamma }\rho \left(\xi \right)dS(\xi )]\\ & =:S_{\Gamma }^w\rho .\end{array}$$

Full text

††journal: J. Comput. Phys.

1 Introduction and problem formulation

A variety of problems of mathematics and physics require the computation of a biharmonic potential subject to Dirichlet boundary conditions. The pure bending problem for an isotropic and homogeneous thin clamped plate is a classical application. Another application is the computation of a $C^{1}$ extension of a given function from its domain of definition to a larger, enclosing domain (we discuss these applications further in section 2.1).

The Dirichlet problem is given as follows. For a domain $D$ with boundary $\Gamma$, find a function $w$ such that

[TABLE]

where $f$ and $g$ are continuous functions defined on $\Gamma$.

The use of standard finite difference methods for the solution of 1, 2 and 3 is complicated greatly by the fact that the differential equation is fourth order. For instance, the resulting linear system for a discretization with $N$ nodes in each dimension would have a condition number proportional to $N^{4}$, which poses several concerns for obtaining high accuracy solutions for large problems.

Integral equation methods, on the other hand, have many advantages for such problems. Because 1, 2 and 3 is homogeneous, the resulting integral equation is defined on the boundary alone and there is a reduction in the dimension of the problem. Complex geometries are handled more easily by an integral equation and, with appropriate choice of representation, the discrete problem tends to be as well conditioned as the underlying physical problem, independent of the system size [1]. One challenge for integral equation methods is that the resulting linear systems are dense. However, there are many well developed fast algorithms for the solution of these systems, most descending from the fast multipole method (FMM) [2].

Integral representations for the solution of 1, 2 and 3 have been developed previously. In particular, the problem is addressed in Peter Farkas’ thesis [3] and the method presented there has been extended to three dimensions in [4]. The integral kernels derived in [3] are taken to be linear combinations of derivatives of the fundamental solution of the biharmonic problem. Assuming the boundary is a smooth curve, the combinations are chosen to maximize the smoothness of the integral kernel as a function on the boundary (for smooth domains). However, the integral kernels derived for 1, 2 and 3 have a leading order singularity of $r^{-2}$ on a domain with a corner. Because of this singularity, designing quadrature rules for discretizing the integral equation is difficult for domains with corners. Furthermore, the resulting discretized system has large condition numbers for domains whose boundaries have high curvature.

For the related problem of two dimensional steady Stokes flow, the stream function formulation results in a biharmonic equation with the gradient of the biharmonic potential specified on the boundary. Let $w$ be the stream function for Stokes flow with no slip boundary conditions, then

[TABLE]

for appropriately chosen functions $f$ and $g$. Over the past century, much work has been done to develop integral representations for the biharmonic problem in this setting, as well as the similar setting of the Airy stress function formulation of the plane theory of elasticity [5, 6, 7, 8]. The representations given in the above references typically have more benign singularities than the representation presented in [3]. In particular, the representation used in this paper, taken from [5, 6], has a leading order singularity of $r^{-1}$ on domains with corners. Moreover, this representation (and others from the above references) can be expressed in terms of Goursat functions, allowing for a convenient representation of the stream function. Because of these advantages, we choose to adapt the representation of [5, 6] to solve 1, 2 and 3.

This adaptation is not immediate. First, in two dimensional Stokes flow, the physical quantities of interest are derivatives of the biharmonic potential $w$ and not $w$ itself; the representation of $w$ from [5, 6] is not necessarily single-valued. Second, in converting the boundary conditions 2 and 3 into the boundary conditions 5 and 6, the data is differentiated along the curve so that the original boundary condition is only met up to a constant. These issues are addressed here, with particular attention paid to the case of multiply connected domains. More precisely, we will show that the desired (and uniquely defined) potential can be expressed in terms of (possibly) multi-valued Goursat functions.

The rest of the paper is organized as follows. In section 2, we present some mathematical preliminaries, including the notation used throughout the paper, a review of the Farkas integral representation, and a review of the completed layer potential representation for solving 4, 5 and 6 in terms of the Goursat functions. In section 3, we explain how to adapt the Stokes layer potentials for the Dirichlet problem, present an integral representation for solving 1, 2 and 3, and prove the invertibility of the resulting integral equation. We outline the numerical tools we used to solve this integral equation and present some numerical results in section 4. In section 5, we provide some concluding remarks and ideas for future research.

2 Preliminaries

The notation for the following concepts can be cumbersome and an attempt has been made to stay consistent. Vector-valued quantities are denoted by bold, lower-case letters (e.g. ${\boldsymbol{h}}$), while tensor-valued quantities are bold and upper-case (e.g. $\mathbf{T}$). Subscript indices of the non-bold character (e.g. $h_{i}$ or $T_{ijk}$) are used to denote the entries within a vector or tensor. We use the standard Einstein summation convention, i.e., there is an implied sum taken over the repeated indices of any term. The vectors ${\boldsymbol{x}}$ and ${\boldsymbol{y}}$ are reserved for spatial variables in ${\mathbb{R}}^{2}$, while $z$ and $\xi$ are reserved for spatial variables in ${\mathbb{C}}$. We make the standard identification between points in ${\mathbb{R}}^{2}$ and points in ${\mathbb{C}}$, i.e. the point ${\boldsymbol{x}}=(x_{1},x_{2})^{\intercal}\in{\mathbb{R}}^{2}$ is equivalent to the point $z=x_{1}+ix_{2}$, and we switch between the two notions implicitly in much of what follows. For integration, the symbol $dS$ is used to denote an integral with respect to arc length and the symbol $d\xi$ is used to denote a complex contour integral. Script letters $\mathcal{X}$, $\mathcal{Y}$, and $\mathcal{Z}$ are reserved for Banach spaces. $I_{\mathcal{X}}:\mathcal{X}\to\mathcal{X}$ denotes the identity operator on $\mathcal{X}$.

Let $D$ will denote a bounded, possibly multiply-connected, domain in ${\mathbb{R}}^{2}$ with a smooth boundary $\Gamma$ (unless otherwise noted). For a domain with $N$ holes, we will denote the outer boundary by $\Gamma_{0}$ and the boundary of each hole by $\Gamma_{1},\ldots,\Gamma_{N}$, so that $\Gamma=\cup_{i=0}^{N}\Gamma_{i}$. Let ${\boldsymbol{n}}({\boldsymbol{x}})$ denote the outward unit normal and $\boldsymbol{\tau}({\boldsymbol{x}})$ the positively-oriented unit tangent for ${\boldsymbol{x}}\in\Gamma$. If we need to distinguish between the exterior and interior of $\Gamma$, we will let $D^{-}=D$ denote the interior and $D^{+}={\mathbb{R}}^{2}\setminus(D\cup\Gamma)$ denote the exterior.

2.1 Applications of the biharmonic Dirichlet problem

Consider the pure bending of an isotropic and homogeneous thin clamped plate. In the Kirchoff-Love theory, the vertical displacement of the plate, $w$, satisfies the equations

[TABLE]

where $D\subset\mathbb{R}^{2}$ represents the midline of the thin plate, $\Gamma$ is its boundary, and $q$ is the transverse load applied to the plate. Using standard techniques, the above problem can be reduced to a homogeneous biharmonic problem of the form 1, 2 and 3.

In a recent paper, [9], it was shown that the solution of polyharmonic Dirichlet problems can be used as part of the solution of inhomogeneous PDEs on complex geometries. We will breifly review this procedure here.

Consider the Poisson equation

[TABLE]

A particular solution, $v$, which satisfies 10 can be obtained from the formula

[TABLE]

where $\Omega$ is some domain such that $D\subset\Omega$ and $\tilde{f}$ is a function defined on $\Omega$ which satisfies $\tilde{f}|_{D}=f$.

There are rapid methods for evaluating the integral 12 in the case that $\Omega$ is a box. However, it is unclear how best to define (and compute) the values of $\tilde{f}$ on $\Omega\setminus D$, particularly such that $\tilde{f}$ is smooth across the boundary of $D$. One approach is to compute the extension as the solution of a homogeneous PDE on the exterior.

Suppose that $w$ solves

[TABLE]

which is a problem of the form 1, 2 and 3. Then, setting $\tilde{f}|_{D}=f$ and $\tilde{f}|_{\Omega\setminus D}=w$ makes $\tilde{f}$ a $C^{1}$ function across $\Gamma$.

In [9], a $C^{0}$ extension was computed as the solution of a Laplace problem on $\Omega\setminus D$. This was found to accelerate the convergence of the Poisson solver over discontinuous extension (i.e. $\tilde{f}$ to be zero outside of $D$). By computing a smoother extension, as in the solution of the problem above, the efficiency and robustness of the Poisson solver could be further improved. For a PDE-based version of this approach, see [10].

2.2 The Farkas integral representation

As mentioned in the introduction, there are existing integral representations for the solution of 1, 2 and 3. In [3], the solution is given as the sum of two layer potentials, i.e.

[TABLE]

where $\sigma_{1}$ and $\sigma_{2}$ are unknown densities.

The integral kernels, $K_{1}^{F}$ and $K_{2}^{F}$ are based on derivatives of the Green’s function for the biharmonic equation. For two points on the plane, ${\boldsymbol{x}}$ and ${\boldsymbol{y}}$, the Green’s function is given by

[TABLE]

Let ${\boldsymbol{r}}={\boldsymbol{y}}-{\boldsymbol{x}}$ and $r=|{\boldsymbol{y}}-{\boldsymbol{x}}|$. Then,

[TABLE]

More explicitly, we have

[TABLE]

On enforcing the Dirichlet boundary conditions for $w$, we obtain the integral equation

[TABLE]

where $\kappa$ denotes the signed curvature as a function on $\Gamma$ and ${\boldsymbol{x}}$ is a point on $\Gamma$. The kernels are given by $K_{11}^{F}=K_{1}^{F}$, $K_{12}^{F}=K_{2}^{F}$,

[TABLE]

For a sufficiently smooth and simply connected domain, the integral equation 20 is invertible. The case of a multiply connected domain is not treated fully in [3], but some of the issues are considered.

As mentioned above, the kernels $K^{F}_{1}$ and $K^{F}_{2}$ are constructed with the goal that the $K^{F}_{ij}$ are as smooth as possible. Suppose that the boundary $\Gamma$ is of class $C^{k}$. Then, the kernels, $K^{F}_{ij}({\boldsymbol{x}},{\boldsymbol{y}})$, are $C^{k-2}$ functions on the boundary for each ${\boldsymbol{y}}\in\Gamma$ [3]. Therefore, on a smooth boundary, these kernels are smooth. However, on a domain with a corner, it is clear from the formula 21 that the kernel $K_{21}^{F}$ has a singularity with strength $r^{-2}$. This singularity, in addition to the term in 20 which explicitly involves the curvature, makes the representation 14 unstable for domains with high curvature (or corners).

2.3 Stokes flow in the plane

The equations of incompressible Stokes flow with no-slip boundary conditions on a domain $D$ with boundary $\Gamma$ are

[TABLE]

where ${\boldsymbol{u}}$ is the velocity of the fluid and $p$ is the pressure. Following standard practice, the velocity ${\boldsymbol{u}}$ can be represented by a stream function $w$. Let

[TABLE]

so that the divergence-free condition, 24, is satisfied automatically. Taking the dot product of $\nabla^{\perp}$ with 23 results in a biharmonic equation for $w$. In particular, $w$ is a biharmonic function which satisfies 4, 5 and 6 with $f=-h_{i}n_{i}$ and $g=h_{i}\tau_{i}$.

2.4 Goursat functions

Goursat showed that any biharmonic function $w$ can be represented by two analytic functions $\phi$ and $\psi$ (called Goursat functions) as

[TABLE]

where $\chi^{{}^{\prime}}=\psi$ and $z=x_{1}+ix_{2}$ [11]. In solving equation 4, 5 and 6, we are interested in $\frac{\partial w}{\partial x_{1}}$ and $\frac{\partial w}{\partial x_{2}}$. Muskhelishvili’s formula [12] gives an expression for these quantites in terms of the Goursat functions as

[TABLE]

We say that a pair of Goursat functions $\phi$ and $\psi$ is equivalent to a Stokes velocity field ${\boldsymbol{u}}$ if the biharmonic function $w$ defined by 27 is such that ${\boldsymbol{u}}=\nabla^{\perp}w$.

The references [5, 6, 7, 8] give many options for the representation of $\phi$ and $\psi$ as layer potentials of a complex density given on the boundary of the domain. Of computational interest, representations for $\phi$ and $\psi$ exist such that enforcing the boundary conditions of 4, 5 and 6 results in an invertible second-kind integral equation (SKIE) for the density.

2.5 Integral representations for Stokes flow in the plane

We will first present the single layer and double layer potentials of Stokes flow in the Stokeslet/stresslet formulation, which may be more familiar. For details concerning these ideas, see, inter alia, [13]. We will then present their equivalent potentials in the classical Goursat function formulation. The reason for doing so is two-fold: first, the Goursat formulation makes it more natural to evaluate the stream function $w$; second, the complex variables-based Goursat formulation is readily adaptable for efficient fast multipole methods.

2.5.1 Stokes layer potentials

The Green’s function $\mathbf{G}$ for the incompressible Stokes equations in free space, or Stokeslet, is given by

[TABLE]

The vector field $u_{i}=G_{ij}\left({\boldsymbol{x}},{\boldsymbol{y}}\right)f_{j}$ represents a Stokes velocity field at ${\boldsymbol{x}}$ due to a point force $\mathbf{f}$ applied at ${\boldsymbol{y}}$. For a continuous distribution of surface forces $\boldsymbol{\mu}$ on a curve $\Gamma$, the induced Stokes field, called a single layer potential, is given by

[TABLE]

The following lemma describes the behavior of the Stokes single layer potential as a function on ${\mathbb{R}}^{2}$, see [13] for details.

Lemma 1**.**

Let $\mathbf{S}_{\Gamma}\mathbf{{\boldsymbol{\mu}}}({\boldsymbol{x}})$ denote a single layer Stokes potential of the form 30. Then, $\mathbf{S}_{\Gamma}{\boldsymbol{\mu}}\left({\boldsymbol{x}}\right)$ satisfies the Stokes equations in ${\mathbb{R}}^{2}\backslash\Gamma$ and $\mathbf{S}_{\Gamma}{\boldsymbol{\mu}}({\boldsymbol{x}})$ is continuous in ${\mathbb{R}}^{2}$.

The Stokes single layer potential has equivalent Goursat functions, $\phi_{S}$ and $\psi_{S}$, which can be expressed in terms of complex layer potentials:

[TABLE]

where $z=x_{1}+ix_{2}$, $\xi=y_{1}+iy_{2}$, and $\rho=\mu_{2}-i\mu_{1}$. The stream function $w_{S}$ corresponding to this Goursat function pair is then

[TABLE]

Note that, the velocity field associated with the stream function $w_{S}$ is given by

[TABLE]

Another quantity of physical interest in Stokes flow is the stress tensor $\boldsymbol{\sigma}$; for a Stokes velocity field ${\boldsymbol{u}}$ and pressure $p$, it is given by

[TABLE]

The stress tensor $\mathbf{T}$, or stresslet, associated with the Green’s function $\mathbf{G}$ is given by

[TABLE]

The vector field $u_{i}=T_{ijk}({\boldsymbol{y}},{\boldsymbol{x}})n_{k}f_{j}$ represents the velocity field resulting from a stresslet with strength oriented in the direction ${\boldsymbol{n}}$ at ${\boldsymbol{y}}$. For a distribution of stresslets $\boldsymbol{\mu}$ on a curve $\Gamma$, the induced Stokes field, called a double layer potential, is given by

[TABLE]

The following lemma describes the behavior of the Stokes double layer potential as a function on ${\mathbb{R}}^{2}$, see [13] for details.

Lemma 2**.**

Let $\Gamma$ be a curve, $D^{+}$ denote the exterior of the curve, $D^{-}$ denote the interior, $\mathbf{D}_{\Gamma}{\boldsymbol{\mu}}\left({\boldsymbol{x}}\right)$ denote a double layer Stokes potential of the form 38, and ${\boldsymbol{x}}_{0}\in\Gamma$. Then, $\mathbf{D}_{\Gamma}{{\boldsymbol{\mu}}}\left({\boldsymbol{x}}\right)$ satisfies the Stokes equations in ${\mathbb{R}}^{2}\backslash\Gamma$ and the jump relations:

[TABLE]

In the above, $\oint$ denotes a Cauchy principal value integral and $\mathbf{D}_{\Gamma}^{PV}{\boldsymbol{\mu}}$ denotes the double layer potential with the integral taken in the principal value sense.

The double layer potential above has equivalent Goursat functions, $\phi_{D}$ and $\psi_{D}$, which can be expressed in terms of complex layer potentials:

[TABLE]

where $z,\xi$, and $\rho$ are as above. The stream function $w_{D}$ corresponding to this Goursat function pair is then

[TABLE]

As before, the velocity field associated with the stream function $w_{D}$ is given by

[TABLE]

2.5.2 The completed double layer representation for Stokes

flow

Using the layer potentials described above, we can represent the solution of $\lx@cref{creftypeplural~refnum}{eq:biharmV1},\lx@cref{refnum}{eq:biharmV2}and~\lx@cref{refnum}{eq:biharmV3}$, or equivalently the system 23, 24, and 25, in terms of a density $\boldsymbol{\mu}$ given on the boundary of the domain. The completed double layer representation [6] for the velocity is

[TABLE]

where $\mathbf{W}_{\Gamma}\boldsymbol{\mu}=\int_{\Gamma}\boldsymbol{\mu}\,dS$, and the representation of an equivalent pair of Goursat functions, giving a stream function $w$, can be inferred from the formulas of the previous subsection. When the no-slip boundary conditions are enforced for this representation, the result is an invertible SKIE for the density $\boldsymbol{\mu}$. The reader can refer to [6] for a detailed discussion of the Fredholm alternative for this representation. We summarize it as

Lemma 3**.**

Let ${\boldsymbol{u}}$ be defined as in 46 and ${\boldsymbol{x}}_{0}\in\Gamma$. Then

[TABLE]

For a sufficiently smooth curve $\Gamma$, the operator $\mathbf{K}_{\Gamma}$ is a compact operator on $\mathcal{X}\times\mathcal{X}$, where $\mathcal{X}$ is $L^{2}(\Gamma)$ or $C^{0,\alpha}(\Gamma)$ for $\alpha\in(0,1)$. Further, the integral equation

[TABLE]

is invertible, even for multiply connected domains.

For the above integral equation, the singularities of the integral kernels which define $\mathbf{K}_{\Gamma}$ are at worst order $r^{-1}$, even for a boundary with a corner.

3 Integral equation derivation

We would like to adapt the completed double layer representation for solutions of Stokes flow 4, 5 and 6 to solve the clamped plate problem 1, 2 and 3. Let $f$ and $g$ be the boundary data as in 1, 2 and 3. By computing tangential derivatives of $f$ on each boundary component, we get the following related Stokes problem:

[TABLE]

There are two main issues to be addressed in using the completed double layer representation in this context. First, the representation is designed for Stokes flow, in which the quantities of interest are derivatives of the potential $\tilde{w}$ and not $\tilde{w}$ itself; the representation for $\tilde{w}$ may not be single-valued. We will establish that, in the context of 49, the stream function is necessarily single-valued. We also discuss some numerical issues related to evaluating the stream function. The second issue to address is that the solution $\tilde{w}$ only satisfies the original boundary condition for the value of $\tilde{w}$ up to a constant on each boundary component. In fact, for multiply connected domains, the completed double layer representation is incomplete for the Dirichlet problem 1, 2 and 3. We present a remedy for this issue and provide some physical intuition.

3.1 Single-valued stream functions

To solve the Dirichlet problem 1, 2 and 3, it is necessary to compute a single-valued biharmonic potential. In the case of a multiply connected domain, there is no guarantee that a single-valued stream function exists for a given velocity field.

Consider the following example. Let $(r,\theta)$ denote standard polar coordinates. It is easy to verify that the velocity field ${\boldsymbol{u}}=\frac{1}{r}\hat{r}$ solves the equations of Stokes flow in an annulus centered at the origin. A stream function for this flow is $w=\theta$, which is not single-valued; indeed, there are no single-valued stream functions which generate this flow.

Let $D$ be a multiply connected domain with boundary $\Gamma=\cup_{i=0}^{N}\Gamma_{i}$, as in the previous section. We note that the gradient of a stream function is determined by the velocity field, i.e.

[TABLE]

Therefore, a velocity field has single-valued stream functions if and only if ${\boldsymbol{u}}^{\perp}$ is conservative. Using standard results from multivariable calculus, we can characterize such flows.

Proposition 1**.**

Suppose that ${\boldsymbol{u}}$ is a divergence-free velocity field which is $C^{1}$ on $D$ and continuous on $D\cup\Gamma$. The field ${\boldsymbol{u}}^{\perp}$ is conservative if and only if

[TABLE]

The equalities 51 constitute $N$ linearly independent constraints on the boundary data because the divergence-free condition 24 implies that $\int_{\Gamma}{\boldsymbol{u}}\cdot{\boldsymbol{n}}\,dS=0$. It turns out that these conditions are satisfied when the Dirichlet problem is recast as a Stokes flow 49, as it is easily verified that

[TABLE]

Thus, any stream function $\tilde{w}$ obtained for the Stokes flow 49 is necessarily single-valued.

3.2 Evaluating the stream function

Given compatible boundary data for the velocity field ${\boldsymbol{u}}$, the completed double layer representation for Stokes flow 46 guarantees the existence of a solution density $\boldsymbol{\mu}$ and a corresponding stream function $\tilde{w}$. The Goursat function formula for $\tilde{w}$, see section 2.5.1, is necessarily single-valued, as explained above, but it is not immediately obvious from the formula that this should be true.

The difficulty in the representation of $\tilde{w}$ comes from the part of the stream function corresponding to the double layer potential 44. The second term in the expression for the double layer potential is

[TABLE]

To compute this term, in a naïve numerical implementation, the question of which is the appropriate branch of the logarithm to use would arise at many steps. To avoid this complication, it is possible instead to compute $v_{1}$, up to a constant, as the harmonic conjugate of the function

[TABLE]

We will use this approach to evalute $v_{1}$ numerically. As a result of the Cauchy-Riemann equations, the harmonic conjugate of $v_{2}$, satisfies the following Neumann problem for the Laplace equation:

[TABLE]

It is possible then to use standard integral equation methods to compute $v_{1}$.

Let $v_{1}=S^{L}_{\Gamma}\sigma$, where $S^{L}_{\Gamma}\sigma$ is the single layer potential for Laplace’s equation, given by

[TABLE]

where $\sigma\in\mathcal{X}=C^{0,\alpha}\left(\Gamma\right)$, for some $\alpha\in(0,1)$, is an unknown density (see [1, 14]). Imposing the Neumann boundary conditions results in the following boundary integral equation for $\sigma$:

[TABLE]

where the operator $K^{L}_{\Gamma}$ is compact, so that the integral equation is second kind. For a derivation of this result, see [1].

It is well known that the operator $\frac{1}{2}I_{\mathcal{X}}+K^{L}_{\Gamma}$ has a one dimensional null space. Thus, we choose to solve the above integral equation subject to the constraint $\int_{\Gamma}\sigma\,dS=0$. Furthermore, it is known that solving the Neumann problem subject to the above constraint is equivalent to solving

[TABLE]

where $W_{\Gamma}\sigma=\int_{\Gamma}\sigma\,dS$.

3.3 Making the representation complete

As mentioned above, the solution $\tilde{w}$ of the auxiliary Stokes problem 49 only satisfies the boundary conditions of the original Dirichlet problem 1, 2 and 3 up to a constant on each boundary component. For a simply connected domain, this constant can be recovered from the fact that adding an arbitrary constant to a stream function does not change the velocity field. Thus, in simply connected domains, there is an equivalence in the solutions of 49 and 1, 2 and 3.

To analyze the case of a multiply connected domain, we first consider radially symmetric solutions on an annulus centered at the origin. Let $w\left(r\right)$ be a radially symmetric biharmonic potential. Then $w(r)$ solves the ordinary differential equation (ODE)

[TABLE]

Four linearly independent solutions of this ODE are $1$, $r^{2}$, $\log\left(r\right)$, and $r^{2}\log\left(r\right)$. For each solution, we can compute the associated velocity field ${\boldsymbol{u}}=\nabla^{\perp}w$. By construction, ${\boldsymbol{u}}$ satisfies the continuity condition 24. For the momentum equation 23 to be satisfied, we need that $\Delta{\boldsymbol{u}}$ is a conservative vector field, which is equivalent to the condition that $\int_{\gamma}\Delta{\boldsymbol{u}}\cdot d{\boldsymbol{\ell}}=0$ for any closed loop $\gamma$ in the annulus. For the first three linearly independent solutions, $\Delta{\boldsymbol{u}}=0$ so that $\Delta{\boldsymbol{u}}$ is trivially a conservative vector field. The fourth solution, on the other hand, has $\Delta{\boldsymbol{u}}=\frac{4}{r}\hat{\theta}$. By considering a curve $\gamma$ encircling the origin, we see that $\Delta{\boldsymbol{u}}$ is not a conservative vector field and that any pressure for the velocity field associated with $r^{2}\log\left(r\right)$ is not single-valued.

The function $\frac{1}{8\pi}r^{2}\log\left(r\right)$ is the Green’s function for the biharmonic equation and is the equivalent of a charge for such problems. The above analysis can be extended to show that any solution of the biharmonic equation with net charge cannot be represented as a Stokes velocity field. In simply connected domains, since $\Delta^{2}w=0$, there can be no net biharmonic charge in the domain. For multiply connected domains with genus $N$, the set of stream functions for Stokes velocity fields misses an $N$ dimensional space of solutions, corresponding to biharmonic charges located in the holes of the domain. Following this reasoning, we obtain a complete representation for biharmonic potentials on multiply connected domains by adding $N$ charges, one per each hole of the domain, to the representation for $w$. The details of this approach, and the proof that it is sufficient, is in the next section.

3.4 The integral representation

Following the discussion in the previous two sections, it is now possible to present an integral representation for the Dirichlet problem of the biharmonic equation based on the completed double layer representation for the Stokes problem. We first fix some notation. Let $D$ be a multiply connected domain, with boundary $\Gamma=\cup_{k=0}^{N}\Gamma_{k}$, as in the previous sections. For each boundary component $\Gamma_{k}$, let $D_{k}$ be its interior and $z_{k}$ be a point in $D_{k}$. Then, let the solution $w$ be represented in terms of layer potentials and biharmonic charges as

[TABLE]

where $\rho$ is an unknown density, the $c_{k}$ are unknown constants, the distance from $z$ to $z_{k}$ is $r_{k}=|z-z_{k}|$, and the operators $S^{w}_{\Gamma}$ and $D^{w}_{\Gamma}$ map complex densities to the stream functions corresponding to single and double layer potentials, as defined in section 2.

Remark 1**.**

As discussed in section 3.2, we will only evaluate the operator $D^{w}_{\Gamma}$ up to a constant in our numerical implementation. However, this does not affect the analysis of this section because of the freedom in choosing $c_{0}$.

As before, we can identify a real, vector-valued density $\boldsymbol{\mu}=(\mu_{1},\mu_{2})^{\intercal}$ with $\rho$ by setting $\mu_{2}({\boldsymbol{x}})-i\mu_{1}({\boldsymbol{x}})=\rho(z)$. Let ${\boldsymbol{u}}=\nabla^{\perp}w$ be the velocity field corresponding to the stream function $w$. Then, in terms of $\boldsymbol{\mu}$, we have

[TABLE]

where $\mathbf{S}_{\Gamma}\boldsymbol{\mu}$ and $\mathbf{D}_{\Gamma}\boldsymbol{\mu}$ are the single and double layer potentials for the density $\boldsymbol{\mu}$, as defined in section 2.

Let $\alpha\in\left(0,1\right)$ and $\mathcal{X}=C^{0,\alpha}\left(\Gamma\right)$. Assume that the boundary data for the Dirichlet problem 1, 2 and 3 satisfies $f\in C^{1,\alpha}\left(\Gamma\right)$ and $g\in\mathcal{X}$, a slightly stronger assumption on the regularity of $f$ than given in the original problem statement. Denote the integrals of $f$ around each boundary component by $b_{k}=\int_{\Gamma_{k}}f\,dS$. To solve equation 1, 2 and 3, we impose the boundary conditions on the gradient of $w$ as in 49 on the above representation for $w$, or, equivalently, the no-slip boundary conditions 25 on the above representation for ${\boldsymbol{u}}$ with

[TABLE]

Under the assumptions on $f$ and $g$, the boundary data ${\boldsymbol{h}}\in\mathcal{X}\times\mathcal{X}$.

Let ${\mathbf{B}}{\boldsymbol{c}}({\boldsymbol{x}})$ denote the part of the velocity field due to the charges, i.e.

[TABLE]

Then, due to lemmas 1 and 2, enforcing the boundary condition ${\boldsymbol{u}}({\boldsymbol{x}})={\boldsymbol{h}}({\boldsymbol{x}})$ for each ${\boldsymbol{x}}\in\Gamma$ results in the following boundary integral equation

[TABLE]

To ensure that the values of $w$ are correct on the boundary, further constraints are needed. We impose $N+1$ additional conditions on the value of $w$

[TABLE]

where the constants $b_{k}$ are as defined above. The integral of $w$ about each component can be written in terms of the unknowns as

[TABLE]

Combining equations 67, 68, and 70, we get the following linear system for the unknowns $\boldsymbol{\mu}$ and ${\boldsymbol{c}}$

[TABLE]

where $\mathbf{D}=\left(D_{0},\ldots D_{N}\right)^{\intercal}$, $\mathbf{F}=\left(F_{0},\ldots F_{N}\right)^{\intercal}$, and ${\mathbf{B}}=\left(b_{0},\ldots b_{N}\right)^{\intercal}$.

Proposition 2**.**

The block system 71 is an invertible Fredholm operator.

Proof.

It is simple to show that the linear system 71 is Fredholm. The block which contains $-1/2\mathbf{I}_{\mathcal{X}\times\mathcal{X}}+\mathbf{K}_{\Gamma}$ is Fredholm due to lemma 3. The off-diagonal blocks, denoted by ${\mathbf{B}}$ and $\mathbf{D}$, are trivially compact because either the domain or range of the operator is finite dimensional. Finally, $\mathbf{F}$ is Fredholm because it is a finite-dimensional linear operator. Therefore, the full system is Fredholm.

Due to the Fredholm alternative, it is only necessary to establish the injectivity of the system 71 to prove that it is invertible. It is clear that if $\boldsymbol{\mu}$ and ${\boldsymbol{c}}$ solve equation 71, then the resulting solution, $w$, given by 62, solves the original Dirichlet problem 1, 2 and 3. By construction, $w$ is biharmonic in $D$. Moreover, $w$ satisfies $\frac{\partial w}{\partial\tau}=\frac{\partial f}{\partial\tau}$ and $\frac{\partial w}{\partial n}=g$ on the whole boundary $\Gamma$ and $\int_{\Gamma_{k}}w=\int_{\Gamma_{k}}f$ for each boundary component $\Gamma_{k}$, so that the boundary conditions are satisfied.

In the case that ${\boldsymbol{h}}\equiv{\boldsymbol{0}}$ and ${\boldsymbol{b}}={\boldsymbol{0}}$, we have that $f=g\equiv 0$ for the Dirichlet problem. By the uniqueness of solutions to 1, 2 and 3, this implies that $w\equiv 0$ in $D$. It is, however, less immediate that $w\equiv 0$ implies that $\boldsymbol{\mu}\equiv{\boldsymbol{0}}$ and ${\boldsymbol{c}}={\boldsymbol{0}}$.

For each $k=1,\ldots,N$, let $\tilde{\Gamma}_{k}\subset D$ be a curve which satisfies $n\left(z_{j},\tilde{\Gamma}_{k}\right)=\delta_{jk}$, where $n\left(z,\gamma\right)$ represents the winding number of the curve $\gamma$ about z. Because ${\boldsymbol{u}}=\nabla^{\perp}w$ and $w\equiv 0$ in $D$, we have

[TABLE]

Let ${\boldsymbol{u}}^{\mu}=\mathbf{S}_{\Gamma}\boldsymbol{\mu}+\mathbf{D}_{\Gamma}\boldsymbol{\mu}={\boldsymbol{u}}-{\mathbf{B}}{\boldsymbol{c}}-\mathbf{W}_{\Gamma}\boldsymbol{\mu}$. We observe that ${\boldsymbol{u}}^{\mu}$ corresponds to a Stokes velocity field in $D$ for any $\boldsymbol{\mu}$. Let $p$ be its associated pressure. Then

[TABLE]

Further, a simple calculation shows that

[TABLE]

for $j=1,\ldots,N$. Combining these equations, we conclude that

[TABLE]

Thus $c_{k}=0$ for $k=1,2,\ldots N$.

The first row of the system 71 then reads

[TABLE]

From the invertibility of $-\frac{1}{2}\mathbf{I}_{\mathcal{X}\times\mathcal{X}}+\mathbf{K}_{\Gamma}$, we conclude that $\boldsymbol{\mu}\equiv{\boldsymbol{0}}$. Because $\boldsymbol{\mu}\equiv{\boldsymbol{0}}$ and $c_{k}=0$ for $k=1,\ldots,N$, we get that $w\equiv c_{0}$. It then follows that $c_{0}=0$ as well, proving the injectivity of the system. ∎

4 Results

We first review the existing numerical tools used to compute solutions of the integral equation 71. To discretize the integral equations, we use the Nyström method. We divide the boundary into panels and represent the unknown density and the boundary data by their values at scaled Gauss-Legendre nodes on each panel. Let $n_{p}$ denote the number of Gauss-Legendre panels. We discretize each panel using 16 scaled Gauss-Legendre nodes. Then $n_{d}=16n_{p}$ is the number of discretization points on the boundary. Let ${\boldsymbol{x}}_{j}$ denote the discretization nodes, $w_{j}$ denote the appropriately scaled Gauss-Legendre quadrature weights for smooth functions, and ${\boldsymbol{\mu}}_{j}$ denote the unknown density at ${\boldsymbol{x}}_{j}$. When forming the linear system, we use scaled unknowns, ${\boldsymbol{\mu}}_{j}\sqrt{w_{j}}$, so that the spectral properties of the discrete system with respect to the $l_{2}$ norm are approximations of the spectral properties of the continuous system as on operator on $L_{2}$ (for more on this point of view, see [15]). The integral kernels in this paper are either smooth or have a weak (logarithmic) singularity. For the smooth kernels in the integral representation, we use standard Gauss-Legendre weights appropriately scaled. For kernels with a logarithmic singularity, we use order $20$ generalized Gaussian quadrature rules [16, 17].

After applying the integral rule, we obtain a linear system for the unknowns. This system is typically well-conditioned, but dense. Let ${\mathbf{A}}$ denote the discretized linear system of size $2n_{d}+N+1$ corresponding to the integral equation (71). Let $\kappa({\mathbf{A}})$ denote the condition number of the discretized matrix ${\mathbf{A}}$. For our applications, the system size was modest and we computed the unknowns $\boldsymbol{\mu}$ and $c_{k}$ using Gaussian elimination. For larger applications, the system is amenable to solution by any of a variety of iterative or fast-direct solvers, which we will not review here.

For the visualizations in this section, we evaluate the layer potentials inside the domain, with some points being very close to the boundary. The value of the layer potential can be difficult to evaluate at such points because of the near-singularity in the integral kernel. We use a sixth order quadrature by expansion method [18, 19] to evaluate these integrals efficiently and accurately.

In this section we consider two test cases. The first example is a convergence study for a simply connected domain to demonstrate the order of convergence for the discretized integral equation. We also compare the condition numbers for the discretized linear systems corresponding to our integral representation and the existing integral representation by Farkas [3] for a family of simply connected domains with increasing curvature. For the second example, we demonstrate the order of convergence and compute the Green’s function for a multiply connected domain. For all examples except the computation of the Green’s function, the boundary data $f$ and $g$ are chosen corresponding to a known solution of the biharmonic equation in $D$ given by

[TABLE]

where $q_{j}$ are uniformly chosen from $[0,1]$. Let $w_{\text{comp}}({\boldsymbol{t}})$ denote the computed solution at targets ${\boldsymbol{t}}$ in the interior of $D$, and let $\varepsilon$ denote an estimate for the error given by

[TABLE]

4.1 Simply connected domain examples

Let $D$ denote the interior of a rounded rectangular bar with length $a=1$, height $b=0.5$, and vertices at $(0,0),(a,0),(a,b),(0,b)$. Following the procedure discussed in [20], the corners are rounded using the Gaussian kernel

[TABLE]

with $h=0.05$. The boundary data $f$ and $g$ are chosen corresponding to a known solution $w$, defined as in 77, with four sources ${\boldsymbol{s}}_{j}$ located at

[TABLE]

with $\delta_{i}$ chosen uniformly from $[-0.05,0.05]$.

The potential $w$ is evaluated at targets ${\boldsymbol{t}}_{j}$ in the interior of $D$,

[TABLE]

with $\delta_{i}$ again chosen uniformly from $[-0.05,0.05]$. A sample geometry with sources ${\boldsymbol{s}}_{j}$ and targets ${\boldsymbol{t}}_{j}$ and the error $\varepsilon$ as a function of $n_{d}$ are shown in fig. 1. The convergence study shows that the error decays like a $20$th order convergent scheme.

The integral equation presented in this paper is significantly better conditioned than the existing integral equation discussed in [3] — particularly when the boundary has regions with large curvature. We plot the condition number $\kappa({\mathbf{A}})$ of the discretized system of integral equations for the representations given by both 71 and 20 as a function of the rounding parameter $h$ for the rounded rectangular bar in fig. 2. The maximum curvature of the boundary is directly proportional to $1/h^{2}$. The condition number $\kappa({\mathbf{A}})$ increases linearly with the maximum curvature for integral equation 20, but is independent of the curvature for the integral equation presented in this paper.

4.2 Multiply connected domain - examples

Let $D$ now denote the interior of a multiply connected domain, where the outer boundary $\Gamma_{0}$ is the boundary of the rounded rectangular bar discussed above with rounding parameter $h=0.05$ and the domain has ten circular obstacles $\Gamma_{i}$ with radii $r_{0}=0.04$ and centers located at ${\boldsymbol{x}}_{i}$,

[TABLE]

We will first perform a convergence study, as above, with a known solution $w$ defined in terms of point sources according to 77. We create ten sources, one located inside each obstacle, whose locations are given by

[TABLE]

where $\delta_{i}$ are chosen uniformly from $[-0.5r_{0},0.5r_{0}]$. The potential is then tested at twelve targets located at

[TABLE]

where $\delta_{i}$ are chosen uniformly from $[-0.5r_{0},0.5r_{0}]$.

We observe $20$th order convergence in the error even for this example. The error as a function of the number of discretization points along with a sample geometry are shown in fig. 3. We also plot the field, and the error in evaluating the potential in the volume using a sixth order quadrature by expansion method in fig. 4. We note that the error observed near the boundary is larger than at the targets used for the convergence study; this is a result of the relatively low order of the quadrature by expansion method and could be improved by increasing the number of points on the boundary.

For the final example, we compute the function $w$ which satisfies the PDE

[TABLE]

where $\delta_{{\boldsymbol{x}}={\boldsymbol{y}}}$ is the two dimensional radially symmetric Dirac delta function centered at ${\boldsymbol{y}}$ and $\{{\boldsymbol{t}}_{j}\}$ are defined in 78. This function describes the vertical displacement of an isotropic and homogeneous thin clamped plate with a transverse load given by point forces at the points ${\boldsymbol{t}}_{j}$. It is also, by definition, a linear combination of the domain Green’s function $G^{D}$, as in

[TABLE]

To compute $w$, we first obtain a particular solution $w_{p}$ which satisfies the PDE in the volume and add to it the solution of a homogeneous problem $w_{h}$ to fix the boundary conditions. We have

[TABLE]

where

[TABLE]

and $w_{h}$ satisfies the following homogeneous biharmonic equation,

[TABLE]

We plot the computed solution in fig. 5.

5 Conclusion

We have presented an integral representation for the biharmonic Dirichlet problem which is stable for domains which have a boundary with high curvature and is applicable to domains which are multiply connected. The representation is based on converting the Dirichlet problem into a problem with velocity boundary conditions, so that classical representations for the velocity boundary value problem can be used. While the technique of [3] — in which integral kernels are chosen by optimizing over the derivatives of an appropriate Green’s function — is general and powerful, the spectral properties of the resulting operator are undesirable for boundaries with high curvature or a corner. Indeed, it seems intuitive that all direct representations for the biharmonic Dirichlet problem should suffer in some way: such an approach asks that one of the integral kernels be singular enough to result in a second kind Fredholm equation for the value of the layer potential and smooth enough to result in a first kind Fredholm equation for the normal derivative of the layer potential.

While some of the above is specific to the biharmonic equation, in particular the use of Goursat functions, it is reasonable to expect the approach to generalize to other high order elliptic problems as well. In particular, there are representations for the modified Stokes equations which are analogues of the completed double layer representation used here [13]. The extension of this method to three dimensions is a topic of ongoing research and will be reported at a later date.

Acknowledgments

M. Rachh’s work was supported by the U.S. Department of Energy under contract DEFG0288ER25053, the Office of the Assistant Secretary of Defense for Research and Engineering and AFOSR under NSSEFF Program Award FA9550-10-1-0180, and the Office of Naval Research under award N00014-14-1-0797/16-1-2123. T. Askham’s work was supported by the U.S. Department of Energy under contract DEFG0288ER25053, the Office of the Assistant Secretary of Defense for Research and Engineering and AFOSR under NSSEFF Program Award FA9550-10-1-0180, and the Office of the Assistant Secretary of Defense for Research and Engineering and AFOSR under award FA95550-15-1-0385. The authors would like to thank Leslie Greengard for suggesting this problem and both Leslie Greengard and Shidong Jiang for many useful discussions.