Modeling mixed boundary conditions in a Hilbert space with the complex variable boundary element method (CVBEM)

The Laplace equation that results from specifying either the normal or tangential force equilibrium equation in terms of the warping functions or its conjugate can be modeled as a complex variable boundary element method or CVBEM mixed boundary problem. The CVBEM is a well-known numerical technique that can provide solutions to potential value problems in two or more dimensions by the use of an approximation function that is derived from the Cauchy Integral in complex analysis. This paper highlights three customizations to the technique.•A least squares approach to modeling the complex-valued approximation function will be compared and analyzed to determine if modeling error on the boundary can be reduced without the need to find and evaluated additional linearly independent complex functions.•The nodal point locations will be moved outside the problem domain.•Contour and streamline plots representing the warping function and its complementary conjugate are generated simultaneously from the complex-valued approximating function.


Method details
The CVBEM has been developed by [1] for the solution of general problems involving Laplace or Poisson equations in multiple dimensions. Modeling begins with a simple closed contour of straight line segments. In a two dimensional complex plane, let V be bounded by a simple closed contour, G, such that By defining (k + 1) equidistant nodal points in each G j such that z j,1 and z j,k+1 are the endpoints of G j , the global nodal coordinates are related to local nodal coordinates by z j,1 = z j and z j,k+1 = z j+1,1 = z j+1 . Fig. 1 shows the global and local nodal numbering conventions. If one defines complex numbers v ji at each node z ji , then degree k complex polynomials N k j ðzÞ are uniquely defined on each boundary element G j .
A global trial function of order k is defined by where d j ¼ 1 z 2 G; 0 otherwise; G k (z) is continuous on G and lim maxjG j j ! 0 G k ðzÞ ¼ vðzÞ: It is assumed that v(z) is analytic on G [ V and that each v ji ¼ vðz ji Þ. Along the boundary G, or exterior to the problem domain union boundary, there are defined n nodal points. For development purposes, the n nodes are assumed defined on G [8]. Later, we will move the nodes outward away from the boundary to demonstrate an addition degree of freedom. The simple closed contour, G, in Fig. 2 is divided into n boundary elements, G jÀ1 , G j ,. . .,G n . For each boundary element, an interpolating polynomial will be used to create a piecewise continuous global interpolation function. In Fig. 2, the boundary, G, is ''severed'' at s = 0 and in the positive direction spans until s = L, the arc length of G. In Fig. 3, the boundary is ''flattened'' and the piecewise function presented. Here, k = 1 is chosen, and the complex polynomials N k j ðzÞ are uniquely defined as first order linear functions. The piecewise function of Fig. 3 is continuous on the boundary G for all z 2 G. The basis function will be used to define a linear global trial function, which is the sum of all nodal basis functions multiplied by a corresponding complex coefficient, w j , the nodal point j value of the function being approximated.
[ ( F i g . _ 1 ) T D $ F I G ] Element End Node Element Interior Node [ ( F i g . _ 2 ) T D $ F I G ] Consider the approximation functionv k ðzÞ defined bŷ v k ðzÞ ¼ 1 2pi From Eq. (2), the global trial function is substituted into the Cauchy Integral formula and integrated over a simply connected two-dimensional complex domain, V, with boundary, G, such that, On each G j , define a local coordinate system by It follows that where N k j ðs j Þ ¼ N k j ðz j ðs j ÞÞ, and g j = (z À z j )/(z j+1 À z j ) for z 2 G. Eq. (9) is solved by factoring (s j À g j ) from N k j ðs j Þ. Let N k j ðs j Þ be of the form where the C j are complex constants in the form (a + bi). Division of N k j ðs j Þ by (s j À g j ) gives where R k j ðzÞ is a complex polynomial of degree k À 1, and Note that d j (z) = |z j À z| and u j+1,j (z) is the central angle between points z j+1 , z j , and z. Fig. 4 shows the special case as z approaches G in the limit.
From Eqs. (6), (7), (9) and (10), summation of the complex boundary element contributions from the m boundary element gives In Eq. (13), it is noted that the N k j ðg j Þ have the form of the assumed shape functions on each g j . Letting node z 1 be on the branch cut of the complex logarithm function ln(z À z) such that z 2 V and z 2 G (see Fig. 5), then (13) can be expanded aŝ where P kÀ1 j is a polynomial of degree (k À 1) defined by and ln(z À z j ) is the principal value of the logarithm function. From the continuity of G k (z), it is seen that at the nodal coordinate z j , and that (z À z j ) is a factor as shown in (15). In (14), the N k m term appears due to the circuit around the branch point of the multiple-valued function ln(z À z j ). Letting [ ( F i g . _ 4 ) T D $ F I G ] From (17), it is seen thatv k ðzÞ is continuous over V and has removable singularities at each boundary element endpoint (nodal coordinate z j , j =1, 2, 3, . . ., m).That is R k (z) and P kÀ1 j are continuous complex polynomials, and lim z ! z j ðz À z j Þlnðz À z j Þ ¼ 0; i:e:v k ðz j Þ ¼ R k ðz j Þ Note that sincev k ðzÞ is analytic in V andv k ðzÞ ¼fðzÞ þ iĉðzÞ wherefðzÞ andĉðzÞ are two dimensional potential and stream functions which satisfy the Laplace equation exactly over V. By forcing the approximation values ofv k ðzÞ to be arbitrarily close (within some e) to the boundarycondition values of v(z) on G, then it is guaranteed by the maximum modulus theorem that the approximation of v(z) is bounded by jvðzÞ ÀvðzÞj e, for all z 2 V.
Because the CVBEM results in a two-dimensional function which is an exact solution to the governing partial differential equation on V, convergence ofvðzÞ to v(z) is then achieved on V [ G by forcing convergence on G. This is shown from (3) and (6) by The global trial function is continuous. Thus,ŵðzÞ is analytic in V, allowingŵðzÞ to be used as an approximation function defined almost everywhere (''ae'') inside V as well as exterior to V. This characteristic separates the CVBEM from other approximations techniques. When solved, the CVBEM approximating integral becomes an approximating function of the form where a j and b j are real constants to be determined, and where a 0 , b 0 , a À1 , and b À1 are also real constants to be determined. The method steps [6] start by using the construct (z À z j ) = R j e iuj , defined at each node j, with location z j , Eq. (19) becomeŝ Using Euler's formula of e iu = (cos u + i sin u), the CVBEM approximation function becomeŝ Further evaluation gives Combining terms from Eqs. (21) and (22), Collecting real and imaginary terms yield We can now separate the approximation function into real,f, and imaginary,ĉ, parts: where the potential functions, or real parts, are given bŷ and where the stream functions, or imaginary parts, are given bŷ Recall that ln j includes the effect of the nodal point logarithmic branch cut rotations.

Matrix formulation and model strategy
After the real and imaginary parts of the CVBEM approximation equation have been developed, the next step is to find the constants, a n , and b n , in thef andĉ functions. The two numerical modeling strategies investigated and subsequently compared for both accuracy and efficiency in this paper are collocation and a least squares approach in a Hilbert space. Each has advantages and disadvantages. Our goal is to highlight which modeling strategy works best for mixed boundary value problems with irregular boundaries.

Collocation
For collocation, the approach is to use the CVBEM by juxtaposing Eq. (19) at each nodal point specified on G (in the limit as z approaches G from inside V). Generally, only one nodal value of either f or c is known at each nodal point [4].
Consider the real portion of the CVBEM with one node for collocation point k. The resulting equation from Eq. (26) is: where p n;k ¼ R n;k ðlnðR n;k Þcosu n;k À u n;k sinu n;k q n;k ¼ R n;k ðlnðR n;k Þsinu n;k À u n;k cosu n;k : (32) Evaluate the above with the five necessary potential collocation points on the problem boundary.
This will result in five linearly independent equations for five collocation points on G, The matrix system to be solved is simply Ax = b, where A is a coefficient matrix and b are the known potential values at each of the collocation points: : Once the system has been configured, substitute the coordinates of the collocation points into the second and third column of the coefficient matrix. It is also necessary to calculate the radius and angle of the collocation points from the singleton node, (see Fig. 6).
The final step is to solve the matrix system. Once these values are known, they can be substituted back into the original equation forfðzÞ. Thef function can now be used to approximate all the potential values within the problem domain. Collocation can also be used in the same way to solve for the streamline equation,ĉðzÞ. The most significant change is that instead of solving for a 0 , one must solve for b 0 .

Least squares
For least squares, the approach is similar except now we will use the evaluation points on the known boundary to create an overdetermined system. From Fig. 1, the values, v j ðzÞ are known for either f j ðzÞ or c j ðzÞ or both. We can now create a vector of measurements along the boundary such that values v 1 . . . v m are available. Using n basis functions (z À z j ) ln(z À z j ), we wish to find the complex coefficients, C j ðzÞ such that the approximating function minimizes the sum of squares such that The need for using Eq. (19) becomes apparent when determining the approximate boundary which is associated with the CVBEM approximator functions,v k ðzÞ.

Approximate boundary
In applying the CVBEM to mixed boundary problems, it is necessary to develop an approximate boundary,Ĝ, upon whichv k ðzÞ satisfies the problem boundary conditions. Engineering problems related to stress and strain are adequate candidates for CVBEM analysis. For stress-free boundary conditions,Ĝ is the collection of points defined bŷ whereŵðzÞ ¼fðzÞ þ iĉðzÞ. Also |z| 2 = x 2 + y 2 where |z| is measured from a selected central point in V. If G coincides with G, then necessarilyvðzÞ ¼ vðzÞ on V [ G. The utility of the approximate boundary concept is in the evaluation of the approximation error. Instead of the analysis of abstract error quantities, the goodness of approximation is determined by visually inspecting the closeness-of-fit betweenĜ and G. In those regions, whereĜ deviates substantially from G, additional evaluation points are placed to reduce the approximation errors from using the selected shape functions.

Analysis and numerical results
As an example of the complex variable boundary element method consider the twisting behavior of a homogeneous, isotropic shaft of an arbitrary, but uniform, cross section that is fixed at one end and subjected to a twisting couple at the other end. If the force and deformation behavior is of interest at some location somewhat removed from either end, then the stress and strain characteristics of the cross section as depicted in Fig. 7 are described by either of the following equations [3]: @ 2 cðx; yÞ @x 2 þ @ 2 cðx; yÞ @y 2 ¼ 0; (38) @ 2 fðx; yÞ @x 2 þ @ 2 fðx; yÞ @y 2 ¼ 0: The quantity c(x, y) is the warping function of the cross-section whereas f(x, y) is the conjugate of c(x, y). If the warping function is known over the cross-section, then the out-of-plane warping displacement and the in-plane shear stresses can be calculated from the expressions v ¼ ucðx; yÞ; t xz ¼ mu @cðx; yÞ @x À y ; t yz ¼ mu @cðx; yÞ @y In the above expressions u is the angle of the twist per unit length, m is the shear modulus, and x, y denote the coordinates of a point located from the center of twist. Furthermore, it should be noted that z represents a coordinate axis and should not be confused with the complex variable z = x + iy. If, on the other hand, the problem is posed in terms of the complementary function f(x, y) then the shear stresses are determined from t xz ¼ mu @fðx; yÞ @y À y While the form of Eqs. (38) and (31) are identical, a solution strategy emerges depending on the manner in which the boundary conditions are specified. If the boundary condition of zero normal stress around the perimeter is posed, then a Neumann boundary condition, i.e. specified normal derivative, best describes the problem. In such a case the nonuniform torsion problem is best posed in terms of the warping function, c(x, y). If on the other hand, the problem is best posed in terms of zero shear around the perimeter, then a Dirichlet boundary condition, i.e. specified functions, best describe the problem. In such a case the problem is best posed in terms of the complementary function, f(x, y). While either solution method is well adapted for solid shafts, it is generally more convenient to operate directly with the warping function, c(x, y), rather than its conjugate, f(x, y), for hollow cross-sections.
The following example is used to analyze the CVBEM with established solutions [7] for shaft crosssections of smooth and sharp corner profiles. Consider the torsion of a solid elliptical cross section with major axis a and minor axis b. The shear-stress-free boundary condition can be expressed in terms of the conjugate function f(x, y) expressed on the boundary as fðx; yÞ ¼ The conjugate function f(x, y) as well as the shear stresses can be shown to be fðx; yÞ ¼   Fig. 8a is due to error. Fig.  8b shows some error, but nearly as much. Fig. 9 shows the relative error for the 6-node model, respectively. The CVBEM relative error along the quarter elliptical section boundary is computed as ½fðx; yÞ À fðx; yÞ=fðx; yÞ wherefðx; yÞ is the CVBEM approximate solution and f(x, y) is the exact solution. The quarter model of Fig. 7 was chosen to take advantage of the problem symmetry and to demonstrate the imposition of f boundary conditions along the exterior curved edge and c along the interior straight edge which is the side extending from the origin to the point (a, 0) and the line extending from the origin to the point (0, b) where a and b are 6.25 and 3.75 respectively. Table 1 summarizes the exact and computed warping function and shear-stress values at points in V using the collocation method. Table 2 demonstrates the exact and warping function calculation using least squares for the same number of basis functions, which are the nodes. The graphical depiction of the CVBEM in Fig. 10 uses computer programs MATLAB, Mathematica, and MATLink to model the mixed boundary problem. Fig. 11 is a time analysis of the efficiency of collocation compared to least squares.

Additional information
The complex variable boundary element method (CVBEM) has been shown to be a mathematically sound approach for modeling two-dimensional potential problems [2]. The foundations of the CVBEM method rests in complex variable theory, namely, the Cauchy integral formula. It tells us that if a  [ ( F i g . _ 9 ) T D $ F I G ]    [ ( F i g . _ 1 0 ) T D $ F I G ] [ ( F i g . _ 1 1 ) T D $ F I G ] Error analysis is an area that makes the CVBEM particularly appealing. It is unique among other numerical methods in that the CVBEM approximating function can be evaluated directly to analyze the error of the approximation. This is because the CVBEM develops an exact representation of the modeling error by the determination of an 'approximate boundary' where the CVBEM approximation exactly satisfies the boundary conditions. That is, the approximate boundary is the locus of points where the CVBEM approximation meets the boundary condition values. This approach is in stark contrast to other numerical methods in their analysis of error. Two sources of error are primary in most numerical methods. The two sources consists of error that results from the approximation in solving the governing equation, and error that result in solving the boundary conditions continuously. Popular numerical methods such as finite elements (FEM) and finite differences (FDM) generate both types of errors in modeling potential problems. Model accuracy is usually estimated by comparing the change in results by increasing the number of nodal points. By this the analyst is seeking convergence by showing that the result is both stable and consistent.