## Abstract

A method is presented for the solution of the biharmonic equation through the use of a combination of polynomial and eigenfunction solutions. Results are presented comparing the method with an exact solution. Excellent agreement is obtained for economical solution representations. Calculations, and comparison with established results, of the displacement of an elastic plate under a uniform applied load are presented. Further, results are presented of the calculation of zero Reynolds number flow in a driven cavity.

## 1. Introduction

In this paper, we present a solution method for the classical biharmonic problem, which consists in finding a continuous function *ϕ*(*x*, *y*), with continuous partial derivatives up to fourth order, which satisfies the homogenous biharmonic equation(1.1)and has prescribed values of the function and its normal derivative on the boundary of the domain, e.g. Meleshko (1998).

The biharmonic equation, besides providing a benchmark problem for various analytical and numerical methods, arises in many practical applications. For example, the bending behaviour of a thin elastic rectangular plate, as might be encountered in ship design and manufacture, or the equilibrium of an elastic rectangle, can be formulated in terms of the two-dimensional biharmonic equation, e.g. Timoshenko & Woinowsky-Krieger (1959). Also, Stokes flow of a viscous fluid in a rectangular cavity under the influence of the motion of the walls, can also be described in terms of the solution of this equation, e.g. Pan & Acrivos (1967), Shankar (1993), Srinivasan (1995), Meleshko (1996) or Shankar & Deshpande (2000). A more recent application of the biharmonic equation has been in the area of geometric and functional design, where it has been used as a mapping to produce efficient mathematical descriptions of surfaces in physical space, e.g. Sevant *et al*. (2000) and Bloor & Wilson (2000).

Interest in solutions of the biharmonic equation and their mathematical properties go back over 130 years, and comprehensive reviews of this work have been given by Meleshko (1998, 2003). In his review article, he concentrates upon the method of superposition in which the solution is described in terms of a sum of separable solutions of the biharmonic equation. In other work, Meleshko (1996, 1998) presents results for Stokes flow in a rectangular cavity in which the solution is based upon the sum of terms consisting of the product of exponential and sinusoidal functions, where the coefficients in the series are determined from the requirement that the prescribed boundary conditions are satisfied. The determination of the coefficients necessitates the solution of an infinite system of linear equations, and Meleshko (1998) describes work that has been done in trying to solve this problem, e.g. Meleshko & Gomilko (1997).

Another approach to the solution of the biharmonic equation described by Meleshko (1998) is to use separable solutions consisting of the product of exponential terms with eigenfunctions of the biharmonic equation that satisfy homogeneous boundary conditions, the so-called Papkovich–Fadle (PF) eigenfunctions. The attractiveness of this approach is the decoupling of the two independent coordinate directions, the final solution being the sum of two independent solutions so obtained. However, it has been noted that this method suffers from poor convergence in the series of PF functions to the imposed boundary conditions (Meleshko 1998). Nevertheless, this method has used been by Joseph & Sturges (1978) and also Shankar (1993) to analyse the zero Reynolds number flow in a driven cavity. Using an eigenfunction approach, Gaydon (1965*a*,*b*) has looked at the specific problem in elasticity of a rectangular plate with a general distribution of normal and shear stress over all four edges. Using physical arguments, he determines a stress function that will produce the resultant stresses over the four edges and the prescribed shears at the corners. The remaining distributed normal and shear stresses on each side are self-equilibrating and amenable to analysis using an eigenfunction expansion. Mathematically, the problem is reduced to solving the biharmonic equation for Airy's stress function subject to specified second derivatives (normal and shear stress) on the boundaries.

In this paper, we revisit the use of biharmonic eigenfunctions, and describe a method which overcomes the difficulties previously encountered. In order to achieve this, the solution is made up of two distinct parts, each part being a solution of the biharmonic equation. The first part is a polynomial solution, formed in such a way as to render the boundary conditions for the second part in a form that is amenable to analysis through the use of PF functions. The second part of the solution is itself then made up of the sum of two independent eigenfunction solutions relating to the two independent variables.

In §2, the method is described in detail for a square domain, although the method is easily extended to rectangular domains. In §3, there is a discussion of how symmetry considerations may be used to reduce the problem. In §4, some results are given which show the accuracy of the method on a test problem where it is possible to compare the solution generated by this approach with an exact solution. It will be seen that a very close approximation to the exact solution can be created, with relatively few eigenfunctions. Section 4 also gives some results for the problem of a uniformly loaded plate, and these are compared with the results of Timoshenko & Woinowsky-Krieger (1959). Further, results for the normal reactions along the edges are also compared with the results of Meleshko *et al.* (2001). Finally, results are also given for the calculation of Stokes flow in a driven cavity.

## 2. Solution method

The system with which we are concerned is that of the solution of the biharmonic equation over *Ω*, the unit square, 0≤*x*≤1, 0≤*y*≤1, in (*x*, *y*) space, subject to appropriate boundary conditions over ∂*Ω*, the boundary of *Ω*. The boundary conditions are such that there are no singularities anywhere on ∂*Ω*. The object of the present exercise is to obtain an approximate solution to the governing partial differential equation, which satisfies the boundary conditions and is expressible in closed form. An approximate numerical solution is, of course, straightforward to obtain, but fails to allow the arbitrary level of refinement so essential for the purposes of either interrogation in the geometrical mapping application or analysis of the physical mechanism.

We require an approximate solution for *ϕ*, where(2.1)over *Ω*, and where *ϕ* is subject to regular boundary conditions over ∂*Ω*. In the case of a periodic solution, *ϕ*(*x*, *y*) may be taken to be a periodic function of *y*, say. Consequently, in this case, it is possible in principle to express *ϕ*(*x*, *y*) as an infinite Fourier series in *y* with the coefficients functions of *x* determined by the boundary conditions. An approximate analytic solution that does not satisfy the boundary conditions exactly can be found using a truncated series. In the more general case, where boundary conditions are imposed along all four sides of the domain, the situation is not so simple, and it is with this situation that this paper is concerned.

The boundary conditions to be imposed on the function and its normal derivatives are taken to be(2.2)where continuity at the edges of the domain is enforced and no singularities introduced. Hence at (0,0),(2.3)with equivalent conditions imposed at the other three corners.

Furthermore, regularity of the boundary conditions indicates that the following extra corner conditions must be satisfied, namely,(2.4)and this ensures that the cross-derivatives at the corners are continuous.

In order to develop an approximate analytic solution, it is appropriate that we proceed in two distinct stages. First, a solution, *ϕ*^{c}(*x*, *y*), to equation (2.1) is found, which satisfies a set of ‘corner’ conditions obtained from equations (2.3) and (2.4). The reason for this is associated with the choice of the eigenfunctions for the problem and the nature of the boundary conditions that can be accommodated by them. Then, we seek a supplementary eigenfunction solution *ϕ*^{p}(*x*, *y*) which satisfies a set of homogeneous boundary conditions as described below. The final solution is composed of the sum of *ϕ*^{c}(*x*, *y*) and *ϕ*^{p}(*x*, *y*).

The corner conditions that must be satisfied at (0,0) by *ϕ*^{c}(*x*, *y*) are(2.5)and similarly at the other three corners.

These conditions ensure that the boundary conditions to be satisfied by the eigenfunction solution, *ϕ*^{p}, are zero and have zero normal derivatives at the corners of the domain. Further, the conditions on the cross-derivatives of *ϕ*^{c} are imposed since the eigenfunction solution cannot accommodate non-zero values of cross-derivatives at the corners. A polynomial solution of equation (2.1) is sought for *ϕ*^{c}, but this is not quite as straightforward as might initially have been anticipated owing to the linear dependency that arises on imposing the conditions to determine the coefficients. It is found that the required degree of the polynomial is six to overcome such difficulties.

The general solution of equation (2.1) is expressible in the form , where *z*=*x*+i*y* and is its complex conjugate. The functions *R*, *S*, *T* and *U* are arbitrary twice differentiable functions of their respective arguments. Thus, the real and imaginary parts of *z*^{n} (*n*≥0), and (*n*≥1) give polynomial solutions of degree *n*, which are solutions of the type in which we are interested. For *n*=0 there is just one independent solution, for *n*=1 there are two and for *n*=2 there are three. For higher values of *n*, there are four independent solutions. Thus, for the general degree six polynomial solution of equation (2.1), there are 22 independent polynomial basis solutions of degrees ranging from 0 to 6, denoted by *P*_{k}(*x*, *y*) (*k*=1, …, 22), say. The form of *ϕ*^{c}(*x*, *y*) is taken to be:(2.6)where the *A*_{k} are constants. In addition to the 16 conditions given by equation (2.5), a further six are required to determine all the coefficients *A*_{k}.

It is certainly the case that the problem could be split into symmetric and anti-symmetric parts with a consequent reduction in the number of unknown coefficients in each part. However, the number is still large and the use of a symbolic manipulator, in this case MAPLE, is desirable. In view of this the complete problem is solved directly although simplifications arising from symmetry considerations are discussed briefly in §3.

Shankar (2004) adopts a similar approach. However, he uses only the values of the dependent variable and its normal derivative at the corners to make the problem homogeneous, and the polynomial is a quartic. This leaves the eigenfunctions incapable of meeting the values of the cross-derivative at the corners.

In principle, it would be possible to specify the remaining six conditions required to find the *A*_{k} arbitrarily, providing none of the existing conditions is violated and the resulting 22 equations are linearly independent. However, it is appropriate to take advantage of this freedom in order to make the function *ϕ*^{c}(*x*, *y*) a closer approximation to the exact solution. Generally speaking, a choice of conditions deduced from the given boundary conditions will achieve this end. However, care must be taken to ensure that linear independence is maintained, and the following supplementary conditions are used:(2.7)

Other choices of the conditions are possible, for example, *ϕ*_{yy}^{c} could have been specified at all of the corners, with *ϕ*_{xx}^{c} specified at only one. Note, however, that symmetry in the choice gives rise to linear dependence of the coefficients.

Equations (2.5) and (2.7) give 22 equations, which determine the constants *A*_{k}.

The function *ϕ*^{p}(*x*, *y*) is now defined by(2.8)

Clearly, *ϕ*^{p}(*x*, *y*) must satisfy equation (2.1), and to obtain the complete solution *ϕ*^{p}(*x*, *y*) must satisfy the modified boundary conditions, obtained from equations (2.2) using equation (2.8). These are,(2.9)

It is important to note that following from the form of *ϕ*^{c}(*x*, *y*) at the corners(2.10)

A solution of equation (2.1) for *ϕ*^{p}(*x*, *y*) subject to the boundary conditions (2.9), can be found in terms of the appropriate eigenfunctions for the problem, namely PF functions (e.g. Joseph & Sturges 1975):(2.11)whose eigenvalues *σ*_{n} satisfy(2.12)

We introduce *ϕ*_{1}(*x*, *y*) satisfying equation (2.1) and the boundary conditions(2.13)together with boundary conditions on *y*=0 and 1, taken from equation (2.9),(2.14)Referring to equation (2.10), the necessity for the corner solution *ϕ*^{c} is now evident, in that the boundary conditions at the ends of the domain are satisfied identically by the PF eigenfunctions. This avoids a difficulty which has sometimes restricted their use in the past, namely the failure to meet the boundary conditions of a series solution constructed from these functions (Meleshko 1998).

The complex eigenvalues *σ*_{n} are well tabulated with their asymptotic form known (Joseph 1977; Joseph & Sturges 1978). It follows, owing to the linearity of the problem, and since *ϕ*_{1} is real, that a solution may be written,(2.15)where *Φ*_{n}(*x*,*σ*_{n}) is the eigenfunction corresponding to the particular eigenvalue, *σ*_{n}, and Re denotes the real part. These eigenfunctions satisfy a biorthogonality condition and this proves very useful when the boundary conditions fall into a certain form, and details of this approach are given by Smith (1952). In the present problem, however, the biorthogonality property offers no distinct advantage as the implementation is, of necessity, in the form of a truncated series approximation and yields a set of simultaneous equations for the coefficients *B*_{n} and *C*_{n}.

Imposing the boundary conditions given in equation (2.14), and using equation (2.15) for *ϕ*_{1}, but with the series truncated to *N* terms, yields the following set of equations for the determination of the constants *B*_{n} and *C*_{n}:(2.16)where the ordering based on *n* is in the ascending order of magnitude of the eigenvalues. A point to note is that the above equations are satisfied identically at *x*=0 and 1 owing to the modification of the boundary conditions by the corner solution *ϕ*^{c}(*x*, *y*).

It is convenient to write each of the above equations (2.16) in the generic form,(2.17)

Following the approach of Shankar (1993), by substituting say, for *i*=2 to *m*−1, into equations (2.17) and forming(2.18)then setting(2.19)where the *α*_{r}, *r*=1, …, 4*N*, are the real and imaginary parts of each of the *B*_{n} and *C*_{n}, one obtains 4*N* equations for the 4*N* unknowns by a standard least squares fit, e.g. Press *et al*. (1990). The quantity *Χ* may be used as a measure of the fit of the truncated eigenfunction expansion to the actual function it represents.

The solution for *ϕ*_{1} thus obtained, satisfies approximately the boundary conditions given by equations (2.13) and (2.14). By interchanging *x* and *y* in the above procedure, an approximate solution *ϕ*_{2} to equation (2.1) satisfying the boundary conditions,(2.20)and(2.21)can be obtained. Hence, it can be seen that the function satisfies equation (2.1) and by reference to equations (2.13), (2.14), (2.20) and (2.21), the general boundary conditions given by equation (2.9).

The solution represented by satisfies equation (2.1) and is approximate only in the sense that the boundary conditions given by equations (2.2) are not exactly satisfied.

With a solution to the general regular problem available, solutions can be found to those problems involving singularities in the boundary conditions, such as the classic two-dimensional Stokes flow in a driven cavity, by introducing a function *ϕ*^{s}(*x*, *y*) to reproduce the local singular behaviour and looking for a solution of the form . An example of this will be presented below.

## 3. Symmetry and rectangular domains

There are certain simplifications which can be made in the cases of symmetric or anti-symmetric problems. These arise, not only from the restricted appropriate choice of the eigenvalues and associated eigenfunctions, but also through the more economical form of the polynomial solution *ϕ*^{c}(*x*, *y*). In both cases, the number of linear algebraic equations to be solved for the unknown coefficients is reduced. As already mentioned, any general problem may be split into a combination of a symmetric part and an anti-symmetric part, although two distinct solutions then need to be combined to obtain the complete solution. To see the position regarding the determination of the corner solution, it is more convenient to think of *Ω* as defined by −1/2≤*x*≤1/2, 1/2≤*y*≤1/2.

For the case symmetric in *x* say, equation (2.6) is simplified by the omission of all those solutions not symmetric in *x*, with the result that the number of unknown coefficients to be determined is reduced to 12. The required conditions are obtained as before bearing in mind the shifted domain and the redundancy owing to symmetry. The situation is algebraically more simple for the case anti-symmetric in *x*. Equation (2.6) is simplified by the omission of all those terms symmetric in *x* with the result that the number of unknown coefficients to be determined is reduced to 10 as might be expected. The required conditions are obtained as before and of course complement those used in the symmetric case.

When the domain *Ω* is rectangular, the eigenvalues for the *x* and *y* directions are no longer equal, but a simple rescaling may be used to obtain general values from those determined by equation (2.12) for the unit interval.

## 4. Results

### (a) Comparison with an exact solution

By way of an example, and in order to illustrate this approach, an exact solution of the problem, specified through equation (2.1) and an appropriate choice of boundary conditions given by equations (2.2) may be used to test the method and assess errors in the approximate solution. It is of course important that the exact form is not expressible exactly through the combination of polynomial and PF functions used in the above analysis. An exact solution chosen, *ϕ*^{e}(*x*, *y*), is obtained through a straightforward application of separation of variables and in this example is given by(4.1)

The specific boundary conditions to be applied in the present exercise corresponding to equations (2.2) are extracted straightforwardly from equation (4.1). From these conditions, the corner values to be used in equation (2.5) can be obtained by direct substitution of *x* and *y*. To complete the specification of the conditions required to obtain *ϕ*^{c}(*x*, *y*), equations (2.5) and (2.7) can be evaluated from equation (4.1) by differentiation and direct substitution as before.

For the results presented three values of *m* in equation (2.18) are used, 16, 31, and 61, and for each the value of the maximum error is calculated as a function of the number of eigenfunctions *N*. The results are presented in figure 1 which shows a plot of log_{10} of the maximum error as a function of the number of eigenfunctions. Naturally, since we are using a least-squares fit in equation (2.19) to determine the coefficients of the eigenfunctions, *N*<*m*. It can be seen that there is a systematic reduction in the value of the error as the number of eigenfunctions *N* used is increased. Furthermore, as the number of points *m* used is increased, the overall error decreases.

Thus, the method presented is capable of producing very close approximations to the exact solution of the biharmonic equation over the unit square. Most importantly, the solution is obtained in closed form in terms of simple analytic functions.

### (b) Practical examples

#### (i) Uniformly loaded plate

We now consider a classical problem in elasticity, which is of great practical importance and has attracted attention over a considerable period of time, e.g. Timoshenko & Woinowsky-Krieger (1959). We refer to problem of a uniform rectangular plate of unit length in the *x*-direction, clamped at its edges and bent by uniform pressure *p* applied to one face. The transverse displacement *w* satisfies the equation(4.2)subject to the boundary conditions that *w* and its normal derivatives are zero at the edges, and where *D* is the flexural rigidity of the plate. The double symmetry in both the *x* and *y* directions is used to simplify the problem as outlined above. Furthermore, the aspect ratio of the plate is taken to be *a*.

This problem has been approached by writing *w* as the sum of two parts, one of which satisfies equation (2.1) and the other acts a particular integral giving rise to the right-hand side equation of (4.2). For example,(4.3)or(4.4)introduced by Boobnoff in 1902 and Hencky in 1913, respectively, as referenced by Meleshko (1997). Alternatively, we can write(4.5)where it can be seen that *ϕ* satisfies equation (2.1) and is subject to the boundary conditions(4.6)

In using the present approach to find the biharmonic function, equation (4.3) represents a special case for this particular set of boundary conditions in that no supplementary corner solution is required. However, to more fully illustrate the current general approach, the second or third option are considered in the present paper. In fact, both cases give rise to virtually identical results for *w*.

An approximate solution of this system is found using the above method with *N*=7 and *m*=31. The results for the displacement at the centre of the plate, measured in units of *p*/*D* are shown in table 1 where the results given by Timoshenko & Woinowsky-Krieger (1959) to about three figure accuracy are also shown. Excellent agreement is obtained, and indeed, bearing in mind the accuracy reported above, where comparison with an exact solution was made, the results for the present problem can reasonably be expected to possess high accuracy.

The local behaviour of the deflection and shear forces near the corner points of the plate are also of interest (Meleshko *et al*. 2001). In order to achieve sufficient spatial resolution to capture this distinctive behaviour, it is necessary to increase the value of *m* to 81 with a corresponding increase in the number of symmetric eigenfunctions *N* to 37. Results showing the normal reactions *V*_{x} along the clamped sides *x*=1/2 and *y*=1 are shown in figures 2 and 3 for a plate with aspect ratio *a*=2. These quantities are related to the deflection *w* by relations given by Meleshko *et al*. (2001). Comparison of these results with the results of Meleshko *et al*. (2001) show that the present method resolves the behaviour of the solution in the vicinity of the corners.

#### (ii) Stokes flow in a driven cavity

An important class of applications involve the solution of the biharmonic equation where the boundary conditions contain singularities. Although there are a variety of ways of solving the biharmonic equation in this situation, an important class of solution techniques involve representing the solution as the combination of a singular part and a regular part, e.g. Srinivasan (1995).

As a further example of the use of the method, we consider the problem of steady two-dimensional Stokes flow in a square container, where the lid of the cavity is driven with uniform velocity parallel to itself. The corners of the cavity are taken to be at (0,0) and (1,0), (1,1) and (0,1).

By seeking a solution for a stream function *Ψ*, the problem can be formulated mathematically as follows:(4.7)subject to the boundary conditions,(4.8)

It can be seen that this problem as it stands is not regular and that there must be singularities in the solution at (0,1) and (1,1), which manifest themselves as points at which the vorticity is infinite. These have been well studied and documented (e.g. Moffatt 1964). Srinivasan (1995) has devised accurate techniques for handling the strongest singularities at the corners whilst causing minimum disturbance to the boundary conditions elsewhere. In particular, by splitting the dependent variable so that the leading singular term valid close to the corners, as derived by Moffatt (1964), is explicitly accounted for. The resulting solution is uniformly valid close to the corners. A similar approach is used in this paper in order to generate a problem which is regular, and so amenable to the present method. We write(4.9)where *Ψ*_{s}(*x*, *y*) satisfies the biharmonic equation and represents the local singularities at the corners, and *Ψ*_{r}(*x*, *y*) also satisfies the biharmonic equation subject to regular boundary conditions and can thus be evaluated using the present method. The boundary conditions for *Ψ*_{r}(*x*, *y*) are obtained from the original conditions on *Ψ* modified by *Ψ*_{s}(*x*, *y*) and result in a form corresponding to equations (2.2)–(2.4), thus making the solution for *Ψ*_{r}(*x*, *y*) amenable to the present analysis. In this approach, unlike Srinivasan's, the original boundary conditions have been significantly modified but this is of no consequence as general arbitrary regular boundary conditions can be handled. Following Moffatt (1964), the solution *Ψ*_{s}(*x*, *y*) is of the form(4.10)where *r*_{1} is the distance from the corner (0,1) *r*_{2} is the distance from the corner (1,1), *θ*_{1} is , *θ*_{2} is , and *f*_{1}(*θ*) is of the form(4.11)with the constants *A*_{1}, *B*_{1}, *C*_{1}, and *D*_{1} determined by the boundary conditions on the walls adjacent to the corner (0,1), and similarly for *f*_{2}(*θ*).

The problem is solved using *N*=73 and *m*=75 and the general flow is shown in figure 4 with a more detailed representation of the flow in the corner shown in figure 5. It can be seen that the solution picks up two of the infinite set of corner vortices, and shows general agreement with the results of Shankar (1993), Srinivasan (1995), Meleshko (1996) and Shankar & Deshpande (2000) on the location and strength of the primary vortex, which according to the present analysis is located at (0.5,0.762) and has *Ψ*=−0.10007. Srinivasan's corresponding results give *Ψ*=−0.100076 at (0.5,0.765), where account has been taken of the different dimensions of the cavity.

## 5. Conclusions

This paper describes an approach to the solution of the biharmonic problem in a rectangular domain, in which a supplementary solution (in this case a polynomial) is used to render the general boundary conditions amenable to the use of eigenfunctions satisfying homogenous conditions at opposite sides of the domain. The complete solution is built-up by combining the sum of the two eigenfunction solutions with the supplementary solution.

Such an approach can find application in a whole range of problems that involve the solution of linear partial differential equations.

## Acknowledgments

The authors would like to thank the referees for their helpful comments.

## Footnotes

- Received August 31, 2005.
- Accepted October 31, 2005.

- © 2006 The Royal Society