An Unsteady Two-Dimensional Complex Variable Boundary Element Method

Author(s)
Bryce D. Wilkins,
Joshua Greenberg,
Brittany Redmond,
Alan Baily,
Nicholas Flowerday,
Adam Kratch,
Theodore V. Hromadka,
Randy Boucher,
Howard D. McInvale,
Steve Horton

Abstract

The Complex Variable Boundary Element Method (CVBEM) procedure is extended to modeling applications of the two-dimensional linear diffusion partial differential equation (PDE) on a rectangular domain. The methodology in this work is suitable for modeling diffusion problems with Dirichlet boundary conditions and an initial condition that is equal on the boundary to the boundary conditions. The underpinning of the modeling approach is to decompose the global initial-boundary value problem into a steady-state component and a transient component. The steady-state component is governed by the Laplace PDE and is modeled using the Complex Variable Boundary Element Method. The transient component is governed by the linear diffusion PDE and is modeled by a linear combination of basis functions that are the products of a two-dimensional Fourier sine series and an exponential function. The global approximation function is the sum of the approximate solutions from the two components. The boundary conditions of the steady-state problem are specified to match the boundary conditions from the global problem so that the CVBEM approximation function satisfies the global boundary conditions. Consequently, the boundary conditions of the transient problem are specified to be continuously zero. The initial condition of the transient component is specified as the difference between the initial condition of the global initial-boundary value problem and the CVBEM approximation of the steady-state solution. Therefore, when the approximate solutions from the two components are summed, the resulting global approximation function approximately satisfies the global initial condition. In this work, it will be demonstrated that the coupled global approximation function satisfies the governing diffusion PDE. Lastly, a procedure for developing streamlines at arbitrary model time is discussed.

The Complex Variable Boundary Element Method (CVBEM) procedure is extended to modeling applications of the two-dimensional linear diffusion partial differential equation (PDE) on a rectangular domain. The methodology in this work is suitable for modeling diffusion problems with Dirichlet boundary conditions and an initial condition that is equal on the boundary to the boundary conditions. The underpinning of the modeling approach is to decompose the global initial-boundary value problem into a steady-state component and a transient component. The steady-state component is governed by the Laplace PDE and is modeled using the Complex Variable Boundary Element Method. The transient component is governed by the linear diffusion PDE and is modeled by a linear combination of basis functions that are the products of a two-dimensional Fourier sine series and an exponential function. The global approximation function is the sum of the approximate solutions from the two components. The boundary conditions of the steady-state problem are specified to match the boundary conditions from the global problem so that the CVBEM approximation function satisfies the global boundary conditions. Consequently, the boundary conditions of the transient problem are specified to be continuously zero. The initial condition of the transient component is specified as the difference between the initial condition of the global initial-boundary value problem and the CVBEM approximation of the steady-state solution. Therefore, when the approximate solutions from the two components are summed, the resulting global approximation function approximately satisfies the global initial condition. In this work, it will be demonstrated that the coupled global approximation function satisfies the governing diffusion PDE. Lastly, a procedure for developing streamlines at arbitrary model time is discussed.

Keywords

Complex Variables, Diffusion Equation, Laplace Equation, Complex Variable Boundary Element Method (CVBEM), Numerical Techniques for Partial Differential Equations

Complex Variables, Diffusion Equation, Laplace Equation, Complex Variable Boundary Element Method (CVBEM), Numerical Techniques for Partial Differential Equations

1. Introduction

In the current work, the Complex Variable Boundary Element Method (CVBEM) is extended to modeling applications of the two-dimensional linear diffusion partial differential equation (PDE), ${u}_{xx}+{u}_{yy}={u}_{t}$ . The proposed solution tech- nique for problems governed by this PDE is based on the standard approach of decomposing the global initial-boundary value problem into two components; namely, a steady-state component and a transient component. The steady-state component is governed by the Laplace PDE, $\Delta {u}_{1}=0$ , and is modeled using the Complex Variable Boundary Element Method. The transient component is governed by the two-dimensional linear diffusion PDE (hereafter referred to as

the diffusion equation), $\Delta {u}_{2}=\frac{\partial {u}_{2}}{\partial t}$ , and is modeled by a linear combination of

basis functions that are the products of a two-dimensional Fourier sine series in the spatial variables x and y and an exponential function in the temporal variable t. The global approximation function is the sum $u={u}_{1}+{u}_{2}$ of the approximate solutions from the two components.

The methodology presented in this work is suitable for use in modeling pro- blems in which the initial condition is equal on the boundary to the boundary conditions. That is, this methodology is intended for modeling problems such that $u\left(x,y,0\right)=f\left(x,y\right)$ on $\Gamma $ , where $f\left(x,y\right)$ represents the boundary con- ditions of the global BVP, and $\Gamma $ is the boundary of the problem domain. When this condition is not satisfied, this methodology can still be used, however, the global approximation function will not satisfy the initial condition on $\Gamma $ .

2. Modeling Approach

The diffusion partial differential equation, given by ${u}_{xx}+{u}_{yy}={u}_{t}$ , can be solved by decomposing the global problem into a steady-state component and a transient component, denoted ${u}_{1}$ and ${u}_{2}$ , respectively. The governing PDE of the steady-state problem is the two-dimensional Laplace equation,

$\frac{{\partial}^{2}{u}_{1}}{\partial {x}^{2}}+\frac{{\partial}^{2}{u}_{1}}{\partial {y}^{2}}=\Delta {u}_{1}=0,$ (1)

and the governing PDE of the transient problem is the diffusion equation

$\frac{{\partial}^{2}{u}_{2}}{\partial {x}^{2}}+\frac{{\partial}^{2}{u}_{2}}{\partial {y}^{2}}=\Delta {u}_{2}=\frac{\partial {u}_{2}}{\partial t}.$ (2)

When numerical techniques are employed to approximate the functions ${u}_{1}$ and ${u}_{2}$ , the steady-state and transient problems are solved approximately and are denoted ${\stackrel{^}{u}}_{1}$ and ${\stackrel{^}{u}}_{2}$ , respectively. The global approximation function, denoted $\stackrel{^}{u}$ , is the sum $\stackrel{^}{u}={\stackrel{^}{u}}_{1}+{\stackrel{^}{u}}_{2}$ .

In this work, the CVBEM outcome is denoted as $\stackrel{^}{\omega}={\stackrel{^}{u}}_{1}+i{\stackrel{^}{v}}_{1}$ , where ${\stackrel{^}{u}}_{1}$ represents the CVBEM approximation of the potential function and ${\stackrel{^}{v}}_{1}$ re- presents the CVBEM approximation of the stream function. The approximate transient solution will be labeled as ${\stackrel{^}{u}}_{2}$ .

2.1. Complex Variable Boundary Element Method Solution of the Steady-State Component

The Laplace equation is an elliptic PDE that models potential problems such as ideal fluid flow, groundwater flow, electrostatic potential, and steady-state heat conduction. There are several numerical techniques that have been successfully employed in solving potential problems such as these including the Finite Element Method [1] and the Method of Fundamental Solutions [2] [3] . Another useful numerical technique for approximating the solution to these problems is the CVBEM, which produces a function that satisfies the strong formulation of the Laplace equation. Consequently, the CVBEM is the topic of numerous papers and books [4] - [9] and has recently been extended to modeling ap- plications of the Laplace equation in three and higher spatial dimensions on problem domains of general geometry, such as a unit sphere on a circular domain. For details regarding the theoretical development of the CVBEM in two spatial dimensions, the reader is referred to [10] [11] [12] . For details regarding the theoretical development of CVBEM type approximation functions in three and higher spatial dimensions, see [13] and [14] .

To approximate a solution to the steady-state problem, the CVBEM is applied to the Dirichlet boundary conditions of the global BVP. Complex polynomials are used in the current work as the family of basis functions in the CVBEM formulation, however, any family (or combination of families) of analytic basis functions could be used. In fact, the basis functions only need to be analytic throughout the problem domain. Collocation of the CVBEM approximation function with the specified global boundary conditions is used to determine the coefficients of the linear combination of the CVBEM approximation function, however, it is noted that other techniques, such as least squares minimization [15] , could also be used to determine the coefficients.

The CVBEM approximation function, $\stackrel{^}{\omega}$ , is a linear combination of the form

$\stackrel{^}{\omega}\left(z\right)={\displaystyle \underset{k=1}{\overset{p}{\sum}}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{c}_{k}{g}_{k}\left(z\right),$ (3)

where
${c}_{k}$ is the k^{th} complex coefficient,
${g}_{k}\left(z\right)$ is the k^{th} member of the family of basis functions being used in the approximation, and p is the number of terms to be used in the linear combination of the approximation function.

The basis functions that are used in the CVBEM approximation function are obtained from complex variable functions that are analytic at least throughout the problem domain. Computational issues may arise depending on the choice of basis functions used in the CVBEM approximation function. For example, basis functions involving branch cuts, such as complex logarithms or reciprocals of complex monomials, among other types of functions, include considerations of positioning branch cuts to lie outside of the problem region in order to enlarge the area of applicability of the CVBEM approximation. Procedures for handling these branch cuts have been well-documented in several papers and books including, [5] and the most recent book [12] and so are not repeated here.

Wherever the basis functions are analytic, it follows from the Cauchy- Riemann equations that both the real and imaginary components of the resulting CVBEM approximation function are harmonic, and consequently, satisfy the two-dimensional Laplace equation. Hence, both the real and the imaginary parts of the CVBEM approximation function, as well as a linear combination of both parts, can be used as a Laplace solver. When entire functions, such as complex polynomials, are used as the basis functions, the real and imaginary parts of the CVBEM approximation function satisfy the Laplace equation throughout the plane. However, when functions with branch cuts are used as the basis functions, the real and imaginary components of the CVBEM approximation function only satisfy the Laplace equation where the basis functions are analytic.

The coefficients of the CVBEM linear combination are complex numbers with both real and imaginary parts, resulting in two degrees of freedom per term used in the CVBEM approximation function. Since the global boundary conditions are Dirichlet, the real part of the CVBEM, which represents the approximation of the potential function, is used to satisfy the global boundary conditions. In this paper, collocation of the real part of the CVBEM approximation function with the specified boundary conditions from the global BVP is the technique used to determine the coefficients of the CVBEM linear combination. Con- sequently, it is necessary to specify 2p boundary conditions in order to uniquely determine the 2p degrees of freedom associated with the CVBEM linear com- bination.

Various techniques can be used to specify the locations of the boundary conditions. Depending on the problem situation, it may be advantageous to use a uniform spacing, or a random spacing. In the general case, the boundary conditions should be reasonably uniformly spaced, especially if the underlying potential function is unknown. An algorithm for locating the position of CVBEM nodes is provided in [16] .

To obtain the necessary real and imaginary components of the CVBEM approximation function observe that by Equation (3), the CVBEM approxi- mation function has the form

$\stackrel{^}{\omega}\left(z\right)={c}_{1}{g}_{1}\left(z\right)+{c}_{2}{g}_{2}\left(z\right)+\cdots +{c}_{p}{g}_{p}\left(z\right).$

Substituting ${c}_{k}={\alpha}_{k}+i{\beta}_{k}$ and ${g}_{k}\left(z\right)={\lambda}_{k}\left(x,y\right)+i{\mu}_{k}\left(x,y\right)$ , it can be shown that the real ( ${\stackrel{^}{u}}_{1}$ ) and imaginary ( ${\stackrel{^}{v}}_{1}$ ) parts of the CVBEM approximation function, which represent the potential and stream functions, respectively, are

$\begin{array}{c}{\stackrel{^}{u}}_{1}\left(x,y\right)={\alpha}_{1}{\lambda}_{1}\left(x,y\right)-{\beta}_{1}{\mu}_{1}\left(x,y\right)+{\alpha}_{2}{\lambda}_{2}\left(x,y\right)-{\beta}_{2}{\mu}_{2}\left(x,y\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\cdots +{\alpha}_{p}{\lambda}_{p}\left(x,y\right)-{\beta}_{p}{\mu}_{p}\left(x,y\right)\\ ={\displaystyle \underset{k=1}{\overset{p}{\sum}}}\left({\alpha}_{k}{\lambda}_{k}\left(x,y\right)-{\beta}_{k}{\mu}_{k}\left(x,y\right)\right)\end{array}$ (4)

and

$\begin{array}{c}{\stackrel{^}{v}}_{1}\left(x,y\right)={\alpha}_{1}{\mu}_{1}\left(x,y\right)+{\beta}_{1}{\lambda}_{1}\left(x,y\right)+{\alpha}_{2}{\mu}_{2}\left(x,y\right)+{\beta}_{2}{\lambda}_{2}\left(x,y\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\cdots +{\alpha}_{p}{\mu}_{p}\left(x,y\right)+{\beta}_{p}{\lambda}_{p}\left(x,y\right)\\ ={\displaystyle \underset{k=1}{\overset{p}{\sum}}}\left({\alpha}_{k}{\mu}_{k}\left(x,y\right)+{\beta}_{k}{\lambda}_{k}\left(x,y\right)\right).\end{array}$ (5)

The methodology is suitable for use with Dirichlet boundary value problems. Thus, boundary conditions from the potential function are specified at the collocation points. Coefficients for the CVBEM approximation function are fitted so that for each of the 2p collocation points, the following relationship holds where
$\left({x}_{q}\mathrm{,}{y}_{q}\right)$ represents the q^{th} collocation point and
${u}_{1}\left({x}_{q}\mathrm{,}{y}_{q}\right)$ is the specified boundary condition from the potential function at that point.

$\begin{array}{c}{\stackrel{^}{u}}_{1}\left({x}_{q},{y}_{q}\right)={\alpha}_{1}{\lambda}_{1}\left({x}_{q},{y}_{q}\right)-{\beta}_{1}{\mu}_{1}\left({x}_{q},{y}_{q}\right)+{\alpha}_{2}{\lambda}_{2}\left({x}_{q},{y}_{q}\right)-{\beta}_{2}{\mu}_{2}\left({x}_{q},{y}_{q}\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\cdots +{\alpha}_{p}{\lambda}_{p}\left({x}_{q},{y}_{q}\right)-{\beta}_{p}{\mu}_{p}\left({x}_{q},{y}_{q}\right)\\ ={u}_{1}\left({x}_{q},{y}_{q}\right)\end{array}$

This implies the matrix equation

$\left[{A}_{1}\right]\left\{c\right\}=\left\{{u}_{1}\right\}\mathrm{,}$ (6)

where $\left\{{u}_{1}\right\}$ is a vector containing the specified potential boundary conditions from the global BVP, $\left[{A}_{1}\right]$ is the matrix obtained from evaluating the CVBEM basis functions at each of the collocation points, and $\left\{c\right\}$ is a vector containing the unknown coefficients ${\alpha}_{1}\mathrm{,}{\beta}_{1}\mathrm{,}\cdots \mathrm{,}{\alpha}_{p}\mathrm{,}{\beta}_{p}$ . Once the coefficients are deter- mined by solving the linear system in (6), they can be substituted back into Equation (4), which is the CVBEM approximation of the steady-state potential function. The resulting function can be used to approximate all of the potential values corresponding to the steady-state solution within the problem domain.

Additionally, the calculated coefficients can be substituted back into Equation (5) and can be used to approximate all of the streamline values of the steady- state solution within the problem domain. Notice that it is possible to ap- proximate all of the streamline values within the problem domain without knowing any streamline boundary conditions. That is, the equation for the stream function is a direct product of the CVBEM due to the orthogonality of the real and imaginary components of the CVBEM approximation function. Accomplishing this with real variable domain techniques such as the Finite Element Method would require post-processing involving an additional numerical scheme. It is noted that in higher dimensional problems, the real and imaginary components are not necessarily orthogonal.

2.2. Solution to the Transient Component

The transient component of the global initial-boundary value problem is modeled by the PDE in Equation (2). The boundary conditions of this problem are Dirichlet and are specified to be continuously zero. The initial condition of this problem is specified as the difference between the initial condition of the global initial-boundary value problem and the CVBEM approximation of the steady-state potential function.

Since the initial condition of the global problem is assumed to be consistent on the boundary with the specified global Dirichlet boundary conditions, the difference between the global initial condition and the CVBEM approximation of the steady-state solution is approximately zero on the boundary. In fact, it is only nonzero due to the error of the CVBEM approximation function in satisfying the global boundary conditions. The approximate transient solution is given in Equation (7) and is a linear combination of basis functions that are the products of a two-dimensional Fourier sine series and an exponential function. It will be shown in Section 2.3 that these basis functions satisfy the diffusion PDE.

${\stackrel{^}{u}}_{2}\left(x,y,t\right)={\displaystyle \underset{i=1}{\overset{m}{\sum}}}{\displaystyle \underset{j=1}{\overset{n}{\sum}}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{a}_{i,j}\mathrm{sin}\left(\frac{i\text{\pi}x}{{L}_{1}}\right)\mathrm{sin}\left(\frac{j\text{\pi}y}{{L}_{2}}\right){\text{e}}^{-{\text{\pi}}^{2}\left(\frac{{i}^{2}}{{L}_{1}^{2}}+\frac{{j}^{2}}{{L}_{2}^{2}}\right)t}$ (7)

In Equation (7), ${\stackrel{^}{u}}_{2}$ is the value of the potential quantity that is associated with the unsteady (transient) component of the problem at a particular location and time, x and y are spatial variables, t is the model time, ${a}_{i\mathrm{,}j}$ is a real coefficient, and ${L}_{1}$ and ${L}_{2}$ are the length and width of the rectangular do- main, respectively.

One collocation point is needed for each term of the series in Equation (7) in order to uniquely determine the coefficients of the linear combination. There- fore, it is necessary to specify the initial condition at mn distinct points within the problem domain. In general, these initial condition collocation points should be located reasonably uniformly spaced throughout the problem domain.

Notice that since sine functions are used it the Fourier series, which are zero whenever $x=0$ , $x={L}_{1}$ , $y=0$ , or $y={L}_{2}$ , Equation (7) is zero continuously along the boundary of the rectangular problem domain. Therefore, the boundary conditions of the transient problem are satisfied by the transient approximation function. However, this result is specifically dependent upon the fact that the problem domain geometry is rectangular.

In order to approximately satisfy the initial condition, we consider the function in Equation (7) when evaluated at $t=0$ . This is

${\stackrel{^}{u}}_{2}\left(x,y,0\right)={\displaystyle \underset{i=1}{\overset{m}{\sum}}}{\displaystyle \underset{j=1}{\overset{n}{\sum}}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{a}_{i,j}\mathrm{sin}\left(\frac{i\text{\pi}x}{{L}_{1}}\right)\mathrm{sin}\left(\frac{j\text{\pi}y}{{L}_{2}}\right).$ (8)

The coefficients
${a}_{i\mathrm{,}j}$ are determined so that for each of the mn collocation points, the following relationship holds where
$\left({x}_{r}\mathrm{,}{y}_{r}\right)$ represents the r^{th} do- main collocation point and
$u\left({x}_{r}\mathrm{,}{y}_{r}\mathrm{,0}\right)-{\stackrel{^}{u}}_{1}\left({x}_{r}\mathrm{,}{y}_{r}\right)$ , represents the difference between the global initial condition and the CVBEM approximation of the steady-state potential function, which is the transient initial condition, at
$\left({x}_{r}\mathrm{,}{y}_{r}\right)$ .

$\begin{array}{c}{\stackrel{^}{u}}_{2}\left({x}_{r},{y}_{r},0\right)={\displaystyle \underset{i=1}{\overset{m}{\sum}}}{\displaystyle \underset{j=1}{\overset{n}{\sum}}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{a}_{i,j}\mathrm{sin}\left(\frac{i\text{\pi}{x}_{r}}{{L}_{1}}\right)\mathrm{sin}\left(\frac{j\text{\pi}{y}_{r}}{{L}_{2}}\right)\\ =u\left({x}_{r},{y}_{r},0\right)-{\stackrel{^}{u}}_{1}\left({x}_{r},{y}_{r}\right)\\ ={u}_{2}\left({x}_{r},{y}_{r},0\right)\end{array}$

This implies the matrix equation

$\left[{A}_{2}\right]\left\{a\right\}=\left\{{u}_{2}\right\}\mathrm{,}$ (9)

where $\left\{{u}_{2}\right\}$ is a vector containing the calculated values of the transient initial condition. Additionally, $\left[{A}_{2}\right]$ is the matrix obtained from evaluating the transient solution basis functions at each of the mn domain collocation points, and $\left\{a\right\}$ is a vector containing the unknown coefficients ${a}_{i\mathrm{,}j}$ . Once the coefficients are determined by solving the linear system in (9), they can be substituted back into Equation (7) and can then be used to approximate all of the potential values corresponding to the transient component within the problem domain. This is the approximation of the transient solution.

2.3. Global Approximation Function

The approximation of the global solution, denoted $\stackrel{^}{u}$ , is achieved by summing the approximate solutions to the steady-state and transient subproblems. The global approximation function is

$\begin{array}{c}\stackrel{^}{u}\left(x,y,t\right)={\displaystyle \underset{i=1}{\overset{m}{\sum}}}{\displaystyle \underset{j=1}{\overset{n}{\sum}}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{a}_{i,j}\mathrm{sin}\left(\frac{i\text{\pi}x}{{L}_{1}}\right)\mathrm{sin}\left(\frac{j\text{\pi}y}{{L}_{2}}\right){\text{e}}^{-{\text{\pi}}^{2}\left(\frac{{i}^{2}}{{L}_{1}^{2}}+\frac{{j}^{2}}{{L}_{2}^{2}}\right)t}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+{\displaystyle \underset{k=1}{\overset{p}{\sum}}}\left({\alpha}_{k}{\lambda}_{k}\left(x,y\right)-{\beta}_{k}{\mu}_{k}\left(x,y\right)\right)\end{array}$ (10)

Due to the decaying exponential function in the approximation of the transient solution, Equation (10) satisfies the intuition that the global approximation function should approach the steady-state approximation as $t\to \infty $ .

To show that Equation (10) satisfies the diffusion PDE, it is necessary to show that ${\stackrel{^}{u}}_{xx}+{\stackrel{^}{u}}_{yy}={\stackrel{^}{u}}_{t}$ . Equivalently, since $\stackrel{^}{u}={\stackrel{^}{u}}_{1}+{\stackrel{^}{u}}_{2}$ , it is necessary to show that

$\begin{array}{l}\Delta \left({\stackrel{^}{u}}_{1}+{\stackrel{^}{u}}_{2}\right)=\frac{\partial}{\partial t}\left({\stackrel{^}{u}}_{1}+{\stackrel{^}{u}}_{2}\right)\\ \Delta {\stackrel{^}{u}}_{1}+\Delta {\stackrel{^}{u}}_{2}=\frac{\partial {\stackrel{^}{u}}_{1}}{\partial t}+\frac{\partial {\stackrel{^}{u}}_{2}}{\partial t}\end{array}$

Since ${\stackrel{^}{u}}_{1}$ is a linear combination of harmonic functions, it follows that

$\Delta {\stackrel{^}{u}}_{1}=0$ . Further, since ${\stackrel{^}{u}}_{1}$ is a function of x and y, it follows that $\frac{\partial {\stackrel{^}{u}}_{1}}{\partial t}=0$ . Therefore, it suffices to show that $\Delta {\stackrel{^}{u}}_{2}=\frac{\partial {\stackrel{^}{u}}_{2}}{\partial t}$ . We shall show that a single term of ${\stackrel{^}{u}}_{2}$ satisfies the PDE $\Delta {\stackrel{^}{u}}_{2}=\frac{\partial {\stackrel{^}{u}}_{2}}{\partial t}$ . It follows that every term of the linear

combination also satisfies the PDE.

$\begin{array}{c}\Delta {\stackrel{^}{u}}_{2}=-{a}_{i,j}{\left(\frac{i\text{\pi}}{{L}_{1}}\right)}^{2}\mathrm{sin}\left(\frac{i\text{\pi}x}{{L}_{1}}\right)\mathrm{sin}\left(\frac{j\text{\pi}y}{{L}_{2}}\right){\text{e}}^{-{\text{\pi}}^{2}\left(\frac{{i}^{2}}{{L}_{1}^{2}}+\frac{{j}^{2}}{{L}_{2}^{2}}\right)t}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-{a}_{i,j}{\left(\frac{j\pi}{{L}_{2}}\right)}^{2}\mathrm{sin}\left(\frac{i\text{\pi}x}{{L}_{1}}\right)\mathrm{sin}\left(\frac{j\text{\pi}y}{{L}_{2}}\right){\text{e}}^{-{\text{\pi}}^{2}\left(\frac{{i}^{2}}{{L}_{1}^{2}}+\frac{{j}^{2}}{{L}_{2}^{2}}\right)t}\\ =-{a}_{i,j}{\text{\pi}}^{2}\left(\frac{{i}^{2}}{{L}_{1}^{2}}+\frac{{j}^{2}}{{L}_{2}^{2}}\right)\mathrm{sin}\left(\frac{i\text{\pi}x}{{L}_{1}}\right)\mathrm{sin}\left(\frac{j\text{\pi}y}{{L}_{2}}\right){\text{e}}^{-{\text{\pi}}^{2}\left(\frac{{i}^{2}}{{L}_{1}^{2}}+\frac{{j}^{2}}{{L}_{2}^{2}}\right)t}\\ =\frac{\partial {\stackrel{^}{u}}_{2}}{\partial t}\end{array}$

Since $\Delta {\stackrel{^}{u}}_{2}=\frac{\partial {\stackrel{^}{u}}_{2}}{\partial t}$ , and since $\Delta {\stackrel{^}{u}}_{1}=\frac{\partial {\stackrel{^}{u}}_{1}}{\partial t}=0$ , it follows that

$\Delta \stackrel{^}{u}=\Delta \left({\stackrel{^}{u}}_{1}+{\stackrel{^}{u}}_{2}\right)=\Delta {\stackrel{^}{u}}_{1}+\Delta {\stackrel{^}{u}}_{2}=0+\Delta {\stackrel{^}{u}}_{2}=\frac{\partial {\stackrel{^}{u}}_{1}}{\partial t}+\frac{\partial {\stackrel{^}{u}}_{2}}{\partial t}=\frac{\partial}{\partial t}\left({\stackrel{^}{u}}_{1}+{\stackrel{^}{u}}_{2}\right)=\frac{\partial \stackrel{^}{u}}{\partial t}.$

Therefore, $\stackrel{^}{u}$ satisfies the diffusion PDE.

3. Example Problem

The envisaged example problem is based upon a two-dimensional rectangular spatial domain with side lengths ${L}_{1}=2$ and ${L}_{2}=1$ . The boundary conditions are Dirichlet, and the initial condition is equal on the boundary to the boundary conditions.

3.1. Problem Formulation

The boundary conditions for the global initial-boundary value problem are given by

$\begin{array}{l}u\left(0,y\right)=2-{\left(2y-1\right)}^{2},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}u\left(2,y\right)=2-{\left(2y-1\right)}^{2},\\ u\left(x,0\right)={\left(x-1\right)}^{2},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}u\left(x,1\right)={\left(x-1\right)}^{2}.\end{array}$ (11)

The initial condition for the global initial-boundary value problem is specified as

$u\left(x,y,0\right)=1+{\left(x-1\right)}^{2}-{\left(2y-1\right)}^{2}.$

Notice that the boundary conditions of the global problem are consistent with the initial condition. That is, on the boundary of the problem domain, the initial condition is equal to the specified boundary conditions.

3.2. Global Initial Condition Results

Figure 1 shows the transient approximation function as it converges to the transient initial condition. Table 1 quantifies the error of the results depicted in Figure 1 by listing the maximum absolute and relative errors of the global approximation function in satisfying the global initial condition. The errors are considered for various numbers of basis functions used in the transient ap- proximation function. In each case, 32 terms are used in the CVBEM ap- proximation of the steady-state solution. The maximum error was assessed by evaluating the global approximation function at 2500 uniformly-spaced points within the problem domain and comparing the approximation outcome with the known initial condition.

All computations were done with the use of the computer program Matlab, but graphical displays were generated with the computer program Mathematica using Matlink to import the computational results from Matlab to Mathematica.

3.3. Steady-State Results

By using analytic complex variable basis functions, the two-dimensional CVBEM develops both potential (real) and stream (imaginary) functions that are harmonic

Figure 1. Approximation of the global initial condition using various numbers of basis functions (specified by n) in the transient approximation function. The maximum absolute and relative errors associated with these graphics are listed in Table 1.

Table 1. Maximum absolute and relative error in the approximation of the initial condition for various numbers of basis functions used in the transient approximation function.

and satisfy the two-dimensional Laplace equation throughout the problem domain, as well as on the problem boundary, and in the exterior of the problem domain except, depending on the choice of basis function used in the CVBEM approximation function, at a finite number of branch points and along a finite number of branch cuts, as discussed earlier. Figure 2 depicts a contour plot of the CVBEM approximation of the steady-state situation. The flow net is gene- rated by evaluating the stream function throughout the problem domain.

Since the analytic solution is a potential function (the solution of the BVP is

Figure 2. Contour plot of CVBEM steady-state approximation, ${\stackrel{^}{u}}_{1}\left(x\mathrm{,}y\right)$ . Collocation points are shown in red, and potential and streamlines are depicted. The streamlines were generated by evaluating the conjugate imaginary component of the CVBEM approximation function.

assumed for discussion purposes to be a potential function, or that the assumed potential function solution of the BVP is the best approximation of the target problem using a potential function, see [12] ), the difference between the CVBEM approximation of the potential function and the analytic solution of the target potential function is itself a potential function. The difference between these two functions is the “error” function, and since it is a potential function, it attains its maximum value on the problem boundary.

Since the error function attains its maximum on the boundary, the error of the CVBEM approximation of the potential function can be assessed by ex- amination of the maximum departure between the boundary values of the CVBEM approximation of the potential function and the boundary values of the analytic potential solution. The absolute modeling error is depicted in Figure 3. The maximum absolute error of the CVBEM approximation function using 32 terms was 0.00513, and the maximum relative error was 0.01350.

There are several techniques for reducing the error of the CVBEM ap- proximation function: 1) Additional basis functions could be used in the CVBEM approximation function, which would require using additional collocation points. Then, the additional collocation points could be located in the areas of high modeling error. Adding additional collocation points where the error is greatest will lower the upper bound of the maximum absolute error; 2) Existing collocation points could be relocated to the areas of high modeling error with no change to the current number of basis functions used in the CVBEM appro- ximation; 3) Different basis functions could be used in the formulation of the CVBEM approximation function.

3.4. Global Solution and Streamline Vector Development

Because the resulting global approximation function is a well-defined function that is continuous and has continuous partial derivatives, the usual vector

Figure 3. Absolute error of CVBEM steady-state approximation on problem boundary. 32 terms were used in the CVBEM approximation function. The values between 0 and 2 on the x-axis correspond to the bottom edge of the rectangular domain, the values between 2 and 3 correspond to the right edge, the values between 3 and 5 correspond to the top edge of the rectangular domain, and the values between 5 and 6 correspond to the left edge.

gradient can be determined throughout the problem domain (as well as in the exterior of the problem domain wherever there are no branch points or branch cuts associated with the CVBEM approximation function). To accomplish development of a vector field, for a specific model time, the global appro- ximation function is evaluated at the selected value of model time, resulting in a spatially-variable function. After developing the corresponding vector gradient function, a vector field depicting orthogonal flow vectors throughout the pro- blem domain can be developed.

Streamlines are developed for various model times in Figure 4. While only demonstrated for a select number of model times, it is noted that streamlines could be developed using this procedure for an arbitrary model time. Due to the normalization used in the statement of the governing flow equations, modeling time advancement utilizes small increments in time as displayed in the figures.

It is also noted that both the global solution as well as the vector gradient outcomes are functions that are defined continuously throughout the problem domain, and that only a sampling of points and outcome values are depicted in Figure 4. Furthermore, these outcome results do not depend on some dis- cretization of points specified throughout the problem domain as part of the computational procedure, such as is usually required in domain methods such as the finite difference method, finite element method, finite volume method, or other similar numerical modeling approaches.

4. Limitations

To ensure that the boundary conditions of the global problem are satisfied, the steady-state solution is fitted to the global boundary conditions. Since the global boundary conditions are satisfied exclusively by the steady-state solution, the boundary conditions for the time-dependent part are homogeneous. Therefore,

Figure 4. This figure depicts the evolution of the global approximation function at various model times (specified by t). The streamline vector trajectories are developed using the usual gradient procedures. The approximations in this figure were created using 64 terms in the transient approximation function and 32 terms in the CVBEM approximation function.

only sine functions are used in the Fourier series approximation of the time- dependent problem. By using the Fourier sine series, the solution to the time- dependent problem is constructed to be zero along the boundary of a re- ctangular domain of dimensions ${L}_{1}\times {L}_{2}$ . However, this is a property of the rectangular geometry of the domain and would not necessarily be true for other geometries. Thus, this technique is currently generally limited to problems with rectangular domain geometries or problems with geometries that are the union of rectangular subdomains. Additionally, this methodology is limited in that it will not satisfy the initial condition on the boundary if $u\left(x\mathrm{,}y\mathrm{,0}\right)\ne f\left(x\mathrm{,}y\right)$ for all $\left(x\mathrm{,}y\right)\in \Gamma $ since the boundary conditions for the transient component would no longer be homogeneous, which means that the two-dimensional Fourier sine series would not satisfy the boundary conditions of the transient problem.

5. Conclusions

By resolving the global initial-boundary value problem into a steady-state com- ponent and a transient component, the CVBEM numerical technique can be applied to modeling applications of the two-dimensional linear diffusion PDE on a rectangular domain with Dirichlet boundary conditions. The steady-state component is governed by the Laplace PDE with boundary conditions specified so as to match the global boundary conditions. This component is modeled with the CVBEM Laplace solver. The transient component is governed by the diffusion PDE with boundary conditions specified so as to be continuously zero and with an initial condition that is specified as the difference between the global initial condition and the CVBEM approximation of the steady-state potential solution. The transient component is modeled by using a linear combination of basis functions that are the product of a two-dimensional Fourier sine series in space and an exponential function in time. The two components are modeled separately, and then the results are summed to yield the global approximation function. The global approximation function was shown to satisfy the strong formulation of the diffusion PDE, which is not true of every numerical method.

Since the global solution is a differentiable function, streamlines can be developed at arbitrary model time by evaluating the gradient of the global approximation function. Further, streamlines associated with the steady-state approximate solution can be developed by evaluating the imaginary component of the CVBEM approximation function since the real and imaginary com- ponents are orthogonal. The results suggest that the methodology results in an approximation function that converges. However, convergence is limited by the ill condition of the matrices in Equations (6) and (9).

It is noted that this technique is currently limited to problems that have rectangular domains due to the reliance of the method on the fact that a two- dimensional Fourier sine series can be designed so as to be zero on the boundary of a rectangular domain. Additionally, this technique is limited to problems where the initial condition on the boundary is equal to the boundary conditions. This is required so that the boundary conditions of the transient component can be modeled as being continuously zero.

Finally, the reader interested in other complex variable techniques for modeling problems governed by the diffusion equation is referred to the papers [17] and [18] . In [17] , the authors present a meshless method for modeling both advection and diffusion problems that is based on the moving least squares technique. The approximation method in [17] is especially important for its ability to model the advection process. In [18] , the authors present a complex variable reproducing kernel particle method for solving two-dimensional variable coefficient advection-diffusion problems. The advantage of this technique is that the shape function of a two-dimensional problem is formed with a one-di- mensional basis function.

Cite this paper

Wilkins, B. , Greenberg, J. , Redmond, B. , Baily, A. , Flowerday, N. , Kratch, A. , Hromadka, T. , Boucher, R. , McInvale, H. and Horton, S. (2017) An Unsteady Two-Dimensional Complex Variable Boundary Element Method.*Applied Mathematics*, **8**, 878-891. doi: 10.4236/am.2017.86069.

Wilkins, B. , Greenberg, J. , Redmond, B. , Baily, A. , Flowerday, N. , Kratch, A. , Hromadka, T. , Boucher, R. , McInvale, H. and Horton, S. (2017) An Unsteady Two-Dimensional Complex Variable Boundary Element Method.

References

[1] Patil, P.V. and Prasad, J.K. (2013) Solution of Laplace Equation Using Finite Element Method. International Journal of Science, Spirituality, Business, and Technology, 2, 40-46.

[2] Reutskiy, S.Yu. (2012) The Method of Approximate Fundamental Solutions (MAFS) for Elliptic Equations of General Type with Variable Coefficients. Engineering Analysis with Boundary Elements, 36, 985-992.

https://doi.org/10.1016/j.enganabound.2012.01.003

[3] Barrero-Gil, A. (2012) The Method of Fundamental Solutions without Fictitious Boundary for Solving Stokes Problems. Computers and Fluids, 62, 86-90.

https://doi.org/10.1016/j.compfluid.2012.03.024

[4] Dean, T. and Hromadka, T. (2009) A Collocation CVBEM Using Program Mathematica. Engineering Analysis with Boundary Elements, 34, 417-422.

[5] Hromadka, T.V. and Whitley, R.J. (1998) Advances in the Complex Variable Boundary Element Method. Springer, New York.

https://doi.org/10.1007/978-1-4471-3611-8

[6] Hromadka, T. and Whitley, R. (1996) A New Formulation for Developing CVBEM Approximation Functions. Engineering Analysis with Boundary Elements, 18, 39-41.

https://doi.org/10.1016/S0955-7997(96)00027-6

[7] Dumir, P. and Kumar, R. (1993) Complex Variable Boundary Element Method for Torsion of Anisotropic Bars. Applied Mathematical Modelling, 17, 80-88.

https://doi.org/10.1016/0307-904X(93)90096-Y

[8] Hromadka, T. and Pardoen, G. (1985) Application of the CVBEM to Non-Uniform St. Venant Torsion. Computer Methods in Applied Mechanics and Engineering, 53, 149-161.

https://doi.org/10.1016/0045-7825(85)90003-9

[9] Hromadka, T.V. and Zillmer, D. (2011) Boundary Element Modeling with Variable Nodal and Collocation Point Locations. Advances in Engineering Software, 43, 96-103.

[10] Hromadka, T. and Guymon, G. (1984) The Complex Variable Boundary Element Method: Applications. International Journal for Numerical Methods in Engineering.

[11] Hromadka, T. and Lai, C. (1987) The Complex Variable Boundary Element Method. Springer-Verlag, New York.

[12] Hromadka, T. and Whitley, R. (2014) Foundations of the Complex Variable Boundary Element Method. Springer, Berlin.

https://doi.org/10.1007/978-3-319-05954-9

[13] Hromadka, T. (2002) A Multi-Dimensional Complex Variable Boundary Element Method. In: Topics in Engineering, Vol. 40, WIT Press, Billerica, MA.

[14] Hromadka, T. and Whitley, R. (2005) Approximating Three-Dimensional Steady-State Potential Flow Problems Using Two-Dimensional Complex Polynomials. Engineering Analysis with Boundary Elements, 29, 190-194.

https://doi.org/10.1016/j.enganabound.2004.07.004

[15] Bohannon, A. and Hromadka, T. (2009) The Complex Polynomial Method with a Least-Squares Fit to Boundary Conditions. Engineering Analysis with Boundary Elements, 33, 1100-1102.

https://doi.org/10.1016/j.enganabound.2009.02.005

[16] Kendall, T., Hromadka, T. and Phillips, D. (2012) An Algorithm for Optimizing CVBEM and BEM Nodal Point Locations. Engineering Analysis with Boundary Elements, 36, 979-984.

https://doi.org/10.1016/j.enganabound.2011.11.008

[17] Wang, J.-F. and Cheng, Y.-M. (2013) New Complex Variable Meshless Method for Advection-Diffusion Problems. Chinese Physics B, 22, Article ID: 030208.

[18] Wang, J.-F. and Cheng, Y.-M. (2013) Analysis of Variable Coefficient Advection-Diffusion Problems via Complex Variable Reproducing Kernel Particle Method. Chinese Physics B, 22, Article ID: 090204.