An Unsteady Two-Dimensional Complex Variable Boundary Element Method

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. How to cite this paper: Wilkins, B.D., Greenberg, J., Redmond, B., Baily, A., Flowerday, N., Kratch, A., Hromadka, T.V., Boucher, R., McInvale, H.D. and Horton, S. (2017) An Unsteady Two-Dimensional Complex Variable Boundary Element Method. Applied Mathematics, 8, 878-891. https://doi.org/10.4236/am.2017.86069 Received: April 5, 2017 Accepted: June 27, 2017 Published: June 30, 2017 Copyright © 2017 by authors and Scientific Research Publishing Inc. This work is licensed under the Creative Commons Attribution International License (CC BY 4.0). http://creativecommons.org/licenses/by/4.0/ Open Access B. D. Wilkins et al.


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 technique 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, ∆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), ∆u 2 = ∂u2 ∂t , and is modeled by a linear combination of basis functions that are the products of a twodimensional 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 problems in which the initial condition is equal on the boundaty to the boundary conditions.That is, this methodology is intended for modeling problems such that u(x, y, 0) = f (x, y) on Γ, where f (x, y) represents the boundary conditions of the global BVP, and Γ 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 Γ.

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, and the governing PDE of the transient problem is the diffusion equation 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 û1 and û2 , respectively.The global approximation function, denoted û, is the sum û = û1 + û2 .In this work, the CVBEM outcome is denoted as ω = û1 + ıv 1 , where û1 represents the CVBEM approximation of the potential function and v1 represents the CVBEM approximation of the stream function.The approximate transient solution will be labeled as û2 .

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 [2] and the Method of Fundamental Solutions [3], [4].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 [5,6,7,8,9,10] and has recently been extended to modeling applications 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 [11], [12], and [13].For details regarding the theoretical development of CVBEM type approximation functions in three and higher spatial dimensions, see [14] and [15].
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 [16], could also be used to determine the coefficients.
The CVBEM approximation function, ω, is a linear combination of the form where c k is the k th complex coefficient, g k (z) 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 at least analytic 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, [6] and the most recent book [13] 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.Consequently, it is necessary to specify 2p boundary conditions in order to uniquely determine the 2p degrees of freedom associated with the CVBEM linear combination.
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 [17].

B.D. Wilkins, et al.
To obtain the necessary real and imaginary components of the CVBEM approximation function observe that by Equation (3), the CVBEM approximation function has the form y), it can be shown that the real (û 1 ) and imaginary (v 1 ) parts of the CVBEM approximation function, which represent the potential and stream functions, respectively, are ( 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 (x q , y q ) represents the q th collocation point and u 1 (x q , y q ) is the specified boundary condition from the potential function at that point.
û1 (x q , y q ) = α 1 λ 1 (x q , y q ) − β 1 µ 1 (x q , y q ) + α 2 λ 2 (x q , y q ) − β 2 µ 2 (x q , y q ) + • • • + α p λ p (x q , y q ) − β p µ p (x q , y q ) = u 1 (x q , y q ) This implies the matrix equation where {u 1 } is a vector containing the specified potential boundary conditions from the global BVP, [A 1 ] is the matrix obtained from evaluating the CVBEM basis functions at each of the collocation points, and {c} is a vector containing the unknown coefficients α 1 , β 1 , . . ., α p , β p .Once the coefficients are determined 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 approximate 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.

Solution to the Transient Component
The transient component of the global initial-boundary value problem is modeled by the PDE in eqn (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 eqn (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.
In eqn (7), û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,j is a real coefficient, and L 1 and L 2 are the length and width of the rectangular domain, respectively.One collocation point is needed for each term of the series in eqn (7) in order to uniquely determine the coefficients of the linear combination.Therefore, 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 , eqn ( 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 eqn (7) when evaluated at t = 0.This is The coefficients a i,j are determined so that for each of the mn collocation points, the following relationship holds where (x r , y r ) represents the r th domain collocation point and u(x r , y r , 0) − û1 (x r , y r ), 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 (x r , y r ).
This implies the matrix equation where {u 2 } is a vector containing the calculated values of the transient initial condition.Additionally, [A 2 ] is the matrix obtained from evaluating the transient solution basis functions at each of the mn domain collocation points, and {a} is a vector containing the unknown coefficients a i,j .Once the coefficients are determined by solving the linear system in (9), they can be substituted back into eqn (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.

Global Approximation Function
The approximation of the global solution, denoted û, is achieved by summing the approximate solutions to the steady-state and transient subproblems.The global approximation function is Due to the decaying exponential function in the approximation of the transient solution, eqn (10) satisfies the intuition that the global approximation function should approach the steady-state approximation as t → ∞.
To show that eqn (10) satisfies the diffusion PDE, it is necessary to show that ûxx +û yy = ût .Equivalently, since û = û1 + û2 , it is necessary to show that ∂t Since û1 is a linear combination of harmonic functions, it follows that ∆û 1 = 0. Further, since û1 is a function of x and y, it follows that ∂ û1 ∂t = 0. Therefore, it suffices to show that ∆û 2 = ∂ û2 ∂t .We shall show that a single term of û2 satisfies the PDE ∆û 2 = ∂ û2 ∂t .It follows that every term of the linear combination also satisfies the PDE.

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.

Problem Formulation
The boundary conditions for the global initial-boundary value problem are given by The maximum absolute and relative errors associated with these graphics are listed in Table 1.
The initial condition for the global initial-boundary value problem is specified as 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.

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 approximation function.In each case, 32 terms are used in the CVBEM approximation of the steady-state solution.The maximum error was assessed by evaluating the global approximation function at 2,500 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.

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 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 generated by evaluating the stream function throughout the problem domain.
Since the analytic solution is a potential function (the solution of the BVP is 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 [13]), 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 examination 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 approximation 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

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 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 approximation 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 problem 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 discretization 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.

Conclusions
By resolving the global initial-boundary value problem into a steady-state component 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 B.D. Wilkins, et al. 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 components 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 eqns ( 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.

)Figure 1 :
Figure1: 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 Table1.

BFigure 2 :
Figure 2: Contour plot of CVBEM steady-state approximation, û1 (x, y).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.

Figure 3 :
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.

Figure 4 :
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.

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.