Question? Leave a message!




Boundary element method Fundamentals and Applications

boundary element method fracture mechanics and boundary element method for beginners boundary element method for general nonlinear differential operators
Dr.NaveenBansal Profile Pic
Dr.NaveenBansal,India,Teacher
Published Date:25-10-2017
Your Website URL(Optional)
Comment
7 The boundary element method 7.1 Integral formulation of boundary-value problems In principle, the boundary element method is just another aspect of the finite element method. However, there is sufficient difference to warrant the new name, which was first coined by Brebbia and Dominguez (1977). The underlying idea is that a boundary-value problem such as that in Section 3.4, involving a partial differential equation of the form (7.1) Lu = f in D subject to the boundary condition (7.2) Bu = g on C, can be transformed to an integral equation   (7.3) vLudA = vf dA D D using any weighting function v. However, there are circumstance in which the  integral LudA may be reduced to an integral over the boundary C by use vD of a reciprocal theorem, for example Green’s theorem for potential-type problems (Green 1828) or Somigliana’s 1886 identity (cited by Becker 1992) for elasticity problems. In the case of a homogeneous equation, f ≡ 0, so that the integral equation is taken only on the boundary, thus reducing the dimensions of the problem. If the equation is non-homogeneous, then there is a domain integral in eqn (7.3) and it would appear that the boundary nature is lost. There are techniques to overcome this problem, for example the dual reciprocity method (Partridge et al. 1992, Wrobel 2002), but we shall not discuss these here. The weighting function used in the boundary element method approach is 2 the fundamental solution (Kythe 1996) associated withL. Suppose thatL =∇ . 2 N.B. Throughout this text, we have used the positive definite operator−∇ . However, since the positive definite nature of the operator is unimportant from a 2 boundary element perspective, we use the standard operator∇ here, in common with other boundary element authors.The boundary element method 245 Also, for convenience, we shall suppose that the boundary condition on C is of the form (7.4) u = g(s)on C , 1 ∂u (7.5) q≡ = h(s)on C , 2 ∂n where C = C + C . 1 2 The fundamental solution, sometimes called the free-space Green’s function, 2 associated with∇ in two dimensions is (see Exercise 7.1) 1 ∗ (7.6) u =− lnR, 2π where R is the position vector of the field point Q (i.e. a general point in D) relative to the source point P (i.e. a point at which the solution is sought); see ∗ Fig. 7.1. We notice that the fundamental solution u has a singularity at P. Define a disc D , of radius , circumference C and centre P (see Fig. 7.1),   and consider what happens as → 0. ∗ We apply the integral equation (7.3) to the region D− D with v = u ,  ∗ since in D− D both u and u are well defined. We shall develop the form for  Poisson’s equation in the first instance:   ∗ 2 ∗ u ∇ udA = u fdA. D−D D−D   Then, using the first form of Green’s theorem, we obtain      ∗ ∂u ∂u 2 ∗ ∗ ∗ u∇ u dA + u − u ds = u fdA ∂n ∂n D−D C+C D−D    D Q R n P 2 Ñ u = 0 D e e C 2 q = h n C u = g 1 Fig. 7.1 Small disc D in the region D. 246 The Finite Element Method ∗ ∗ 2 ∗ or, using the notation q≡ ∂u/∂n and q ≡ ∂u /∂n and noting that ∇ u = 0in D− D ,    ∗ ∗ ∗ (7.7) (u q− uq ) ds = u fdA. C+C D−D   Now, on C , ∂/∂n≡−∂/∂R and R = . Also, on C , u(s)= u +   P η (s)and q(s)= q + η (s), where max(η ,η )→0as → 0. 1 P 2 1 2 We write  ∗ I = u qds 1 C  and   ∗ ∗ uq ds = u q ds + I , p 2 C C   where  ∗ I = η q ds. 2 2 C  Then, on C ,         1   I = − lnR (q + η ) ds 1 P 2   2π C  1 ≤ ln(q +η )2π P 2 2π →0as → 0. Also,         ∂ 1   I = η − − lnR ds 2 2   ∂R 2π C  1 1 ≤η 2π 2 2π  →0as → 0. Finally,  1 1 ∗ u q ds = u 2π P P 2π  C  = u . P Hence, as → 0,  ∗ ∗ (u q− uq ) ds→ u P C The boundary element method 247 and we have the integral equation   ∗ ∗ ∗ (7.8) u = (u q− q u)+ u fdA P C D for values of u inside D in termsofvaluesof u and q on the boundary. In Exercise 7.3, we obtain the corresponding integral equations for points on the boundary and for points outside the boundary in the form   ∗ ∗ ∗ (7.9) c u = (u q− q u)+ u fdA, P P C D where ⎧ 1,P ∈ D, ⎨ (7.10) c = α /2π, P ∈ C, P P ⎩ 0,P ∈ D∪ C, and where α is the internal angle between the left- and right-hand tangents at P 1 P. If the boundary is smooth at P,then c = . In the terminology of integral P 2 ∗ ∗ equations, the functions u and q in eqn (7.9) are called kernels of the integral equation (Manzhirov and Polyanin 2008). A result which provides a useful check in the boundary element method is obtained by choosing v = 1 in the first form of Green’s theorem to obtain   ∂u 2 ∇ udA = ds, ∂n D C so that if u satisfies Laplace’s equation, then  (7.11) qds=0. C 7.2 Boundary element idealization for Laplace’s equation We proceed in a manner analogous to that for the finite element method in Chapter 3. We shall approximate the boundary C by a polygon, C ,and it n is the polygon edges which are the boundary elements. We also choose a set of nodes, at which we seek approximations U and Q to u and q, respectively. We shall consider in the first instance the so-called constant element. In such an element, the geometry is a straight line segment with just one node at the centre; see Fig. 7.2. N.B. (i) In this element, the approximation of the function is of a lower order than that for the geometry, and the element is known as a superparametric element. Similarly, we can define subparametric elements (cf. the isoparametric elements of Chapter 4). (ii) There is no requirement for interelement continuity in the boundary element, non-conforming elements are frequently used.248 The Finite Element Method l j j j Fig. 7.2 Constant element. We choose a set of basis functions w (s): j =1,2,...,n to approximate j u and q as follows: n n (7.12) u ˜ = w (s)U and q ˜ = w (s)Q , j j j j j=1 j=1 where s is a measure of the arc length around C (cf. eqn (3.16) for a finite n element approximation). The weighting functions for the constant element are given by  1,j∈ e, (7.13) w (s)= j 0,j∈/ e. We now use point collocation with the approximations (7.12) in eqn (7.9), remembering that for Laplace’s equation f ≡ 0, the curve C is replaced by C n and the boundary point P is chosen to be, successively, the nodes 1,2,...,n: ⎛ ⎞     n n 1 1 ∂ ⎝ ⎠ c U = w (s)Q (−lnR ) ds− (w (s)U ) − (lnR ) ds i i j i i j j i 2π 2π ∂n C C n n j=1 j=1 or 1 2 1 2   n n ∂ α U = (ln R ) ds U − lnR ds Q , i i i j i j ∂n j j j=1 j=1 j =1,2,...,n, where R =R, R (s) is the position vector of a boundary point s relative to i i node i,and α =2πc . i i We can rewrite these equations as n n (7.14) H U + G Q =0,j =1,2,...,n, ij j ij j j=1 j=1 where  ∂ (7.15) H = (lnR ) ds− α δ ij ij i ij ∂n jThe boundary element method 249 j R ij Target element i Base node Fig. 7.3 Target element j relative to base node i. and  (7.16) G =− lnR ds, ij ij j with R =R and where R (s) is the position vector of a point in the target ij ij ij element jrelativetothe base node i; see Fig. 7.3. Finally, then, the system of equations (7.14) may be written in matrix form as (7.17) HU + GQ = 0, where U and Q are vectors of the approximations to the boundary values of u and q, respectively. For properly posed problems, we have boundary conditions of the form of eqns (7.4) and (7.5) so that at any point we know only the value of either u or q, and we partition the matrices in eqn (7.13) appropriately. Suppose that U and 1 Q are vectors of the known values of u and q and that U and Q are vectors 2 2 1 of unknown values. Then we may write eqn (7.17) as     U Q 1 1 1 2 1 2 (7.18) H H + G G = 0, U Q 2 2 and hence we can rearrange it in the form (7.19) Ax = b, where   U 1 2 1 1 2 (7.20) A=H G , b =− H G Q 2 and the unknowns are given by the vector  U 2 (7.21) x = . Q 1250 The Finite Element Method Equation (7.19) is an n× n system of algebraic equations for the n unknowns x . We note here that the system matrix A is densely populated, non-symmetric i and non-positive definite. This is in stark contrast to the stiffness matrix for the finite element method in Section 3.6. The solution of the system (7.19) yields the nodal vectors U and Q,and values at k internal points, k =1,2,...,m, may be obtained using the discrete form of eqn (7.9) for k inside the boundary: ⎡ ⎛ ⎞ ⎛ ⎞ ⎤  n n 1 ∂ ⎣ ⎝ ⎠ ⎝ ⎠ ⎦ U = w (s)U lnR − w (s)Q lnR ds, k j j k j j k 2π ∂n C n j=1 j=1 k =1,2,...,m. We write this in matrix form as ˘ ˘ (7.22) U = HU + GQ, I where   1 ∂ 1 ˘ ˘ (7.23) H = (lnR ) ds and G =− lnR ds. kj kj kj kj 2π ∂n 2π j j The integrals required to evaluate H and G , eqns (7.15) and (7.16), involve ij ij the kernels (∂/∂n)(ln R)and lnR and may be considered to be of two types: 1. Base node not in target element. In this case, all the integrals are non-singular and standard Gauss quadrature can be used. It is also possible to obtain analytic values; see Exercise 7.5. 2. Base node in target element, i = j singular integral. In this case, the integrals are singular and we cannot use a standard Gauss quadrature. We proceed as follows. To evaluate the integrals for G , which have a logarithmic singularity ij (eqn (7.16)), we can use a special logarithmic Gauss quadrature; see Appendix D. However, in this case we can obtain analytic values for the singular integrals, and we shall see how to do that in Section 7.3. The singularity in the integral for H , eqn (7.15), is O(1/R), and such ij integrals usually require special treatment. The term H also requires com- ij putation of the parameter α . Both these difficulties may be overcome in the i case of the Laplace operator. We notice that we can apply our approach to the problem whose unique solution is u≡ 1. The equivalent boundary-value problem is 2 ∇ u=0 in D, u=1 on C, and clearly q =0 on C.The boundary element method 251 T T Hence U=11... 1 and Q=00... 0 , satisfying eqn (7.15) so that HU = 0. Consequently, it follows that n  (7.24) H =− H,i=1,2,...,n, ii ij j=1 and the diagonal terms are found from the sum of the off-diagonal terms, a very convenient result. ˘ ˘ Finally, all integrals in the computation of H and G (eqn (7.21)) are non- jk jk singular and a standard Gauss quadrature can be used. Once all elements of the matrices H and G are evaluated, the matrices A and b are assembled (eqn (7.20)), and the system (7.19) is solved by a suitable routine. This leads to values for U and Q on the boundary, and internal values are obtained using eqn (7.22). We note here that even though the evaluation of the internal values uses the approximate boundary values, the internal values are often more accurate than the computed values on the boundary. However, care is needed for internal points close to the boundary, since the values of R in kj eqns (7.23) can become very small compared with the element length. A useful rule of thumb is that internal points should be no closer to the boundary than one quarter of the minimum length of a boundary element. 7.3 A constant boundary element for Laplace’s equation The element is shown in Fig. 7.3, and in Fig. 7.4 we define some of the geometry of element j, whose length is l . j 2 n j q ij j j 1 R ij d ij i Fig. 7.4 Geometry for the constant element j.252 The Finite Element Method We calculate the coefficients H and G . Suppose that the base node i is ij ij not in the target element j. Then ∂ (lnR )=(gradR· n) ij ij ∂n 1 = R · n ij j R ij cos θ ij = R ij d ij = 2 R ij so that, using Gauss quadrature, we obtain  d ij H = ds ij 2 R j ij  1 l d 1 j ij = dξ 2 2 R (ξ) −1 ij G l d 1 j ij (7.25) ≈ w . g 2 2 R (ξ ) g ij g=1 We note here that if node i is in element j, then n is perpendicular to R and j ij R · n =0, so that H =−α . We shall use eqn (7.24) to compute H and, ij j ii i ii where appropriate, use this result as a check:  G =− lnR ds ij ij j  1 l j =− lnR (ξ) dξ ij 2 −1 G l j (7.26) ≈− lnR (ξ )w . ij g g 2 g=1 Now, to obtain G we use the notation of Fig. 7.5: ii 2 R ii i 1 i Fig. 7.5 Base node in target element.The boundary element method 253  G =− lnR ds ii ii j    1 l l i i =−2 ln η dη 2 2 0  1     1 =−l ln l +lnη dη i i 2 0  1   1 (7.27) =−l ln l − l lnηdη. i i i 2 0    1 In Exercise 7.4, we show that G = l 1− ln l . The coefficients H and ii i i kj 2 G are obtained using Gauss quadrature. In Exercise 7.5, we obtain analytic kj expressions for H and G . ij ij Example 7.1 Solve Laplace’s equation in the unit square subject to the bound- ary conditions u(x,0) = u(x,1) = 1 + x, q(0,y)=−1,q(1,y)=1, with four constant boundary elements as shown in Fig. 7.6. We shall use four-point Gauss quadrature for the integrals, and with the notation of Fig. 7.6 we find R (η ): 0.504798 0.599088 0.835995 1.056389 12 g so that 4 1 (ξ )w =4.426960. g g 2 R 12 g=1 3 42 1 Fig. 7.6 Boundary nodes for the problem in Example 7.1.254 The Finite Element Method Hence, using eqn (7.25),   1 1× 2 H = × 4.426960 = 1.106740. 12 2 By symmetry, H = H . Similarly, 14 12 H =0.927812. 13 Now, H =−(H + H + H ) 11 12 13 14 =−3.140762. Now α = π, so we see that H ≈−α as expected. 1 11 1 The values H (i=1,2,3,4;j =1,2,3,4) may be obtained by symmetry. ij Also, 4 lnR (η )w =−0.669656. 12 g g g=1 Hence, using eqn (7.26),   1 G = − × (−0.669656) = 0.334828. ij 2 By symmetry, G = G . Similarly, 14 12 G =0.038869. 13 Using the result of Exercise 7.3, we have G =1.693147. 11 Hence our overall system matrices are 1 2 3 4 ⎡ ⎤ −3.14076 1.10674 0.92781 1.10674 1 ⎢ ⎥ 1.10674 −3.14076 1.10674 0.92781 2 ⎢ ⎥ H = , ⎣ ⎦ 0.92781 1.10674−3.14076 1.10674 3 1.10674 0.92781 1.10674 −3.14076 4 1 2 3 4 ⎡ ⎤ 1.69315 0.33483 0.03887 0.33483 1 ⎢ ⎥ 0.33483 1.69315 0.33483 0.03887 2 ⎢ ⎥ G = . ⎣ ⎦ 0.03887 0.33483 1.69315 0.33483 3 0.33483 0.03887 0.33483 1.69315 4 Now we apply the boundary conditions and partition H and G accordingly. At nodes 2 and 4 we specify u, and at nodes 1 and 3 we specify q; hence we obtainThe boundary element method 255 2 4 1 3 ⎡ ⎤ 1.10674 1.10674 1.69315 0.03887 ⎢ ⎥ −3.14076 0.92781 0.33483 0.33483 ⎢ ⎥ A = , ⎣ ⎦ 1.10674 1.10674 0.03887 1.69315 0.92781 −3.14076 0.33483 0.33483 3 2 4 1  T b = , 3.31943 −4.97450 3.31943 −1.66594 and solving Ax = b yields 3 2 4 1  T x = 1.907 1.094 −0.001 −0.001 . Hence we have the nodal values 4 1 2 3  T U = 1.5 1.907 1.5 1.094 , 4 1 2 3  T Q = −0.001 1 −0.001 −1 . The exact solution is u=1+ x. We see that the approximate values given above compare very well, given the very coarse approximation. The approximate value 4 of qds is−0.002, which again is close the the exact value, zero (eqn (7.11)). C Equation (7.23) yields the following coefficients for the computation of internal values (cf. eqns (7.25) and (7.26)): G G l d 1 l j kj j ˘ ˘ H ≈ w , G ≈− lnR (ξ )w . kj g kj kj g g 2 4π R (ξ ) 4π g kj g=1 g=1   1 1 We shall obtain the solution at , . 2 2   1 1 At , , 2 2 R (ξ ): 0.659840 0.528107 0.528107 0.659840. 11 g Hence ˘ ˘ ˘ ˘ H =0.24965 = H = H = H 11 12 13 14 and ˘ H =0.999, ij very close to 1; see Exercise 7.7. Also, ˘ ˘ ˘ ˘ G =0.08928 = G = G = G . 11 12 13 14256 The Finite Element Method Now,   4 4 1 1 ˘ ˘ u , = H U + G Q 1j j 1j j 2 2 j=1 j=1 =1.498. N.B. As mentioned earlier, this internal potential value is more accurate than the potential values on the boundary.   1 3 The solution at , is given in Exercise 7.6. 2 4 7.4 A linear element for Laplace’s equation We shall consider an element with two nodes, just as we did for the linear finite element in Section 3.5, illustrated in Fig. 3.7. The geometry and parameters for the element are the same as those for the constant element, see Fig. 7.4. ˘ ˘ In this case we find that the integrals for H , G , H and G are ij ij kj kj taken over the two elements containing node j. The details are obtained in Exercise 7.9. Example 7.2 Consider the problem of Example 7.1 using the linear element of Exercise 7.9. We shall assume that we have a Dirichlet problem with u known at all four nodes. 1 2 3 4 ⎡ ⎤ −1.57081 0.43883 0.69312 0.43883 ⎢ ⎥ 0.43883 −1.57081 0.43883 0.69312 ⎢ ⎥ H =⎢ ⎥ , ⎣ 0.69312 0.43883 −1.57081 0.43883⎦ 0.43883 0.69312 0.43883 −1.57081 1 2 3 4 ⎡ ⎤ 1.50000 0.21460 −0.19315 0.21460 ⎢ ⎥ 0.21460 1.50000 0.21460 −0.19315 ⎢ ⎥ G = . ⎢ ⎥ ⎣ ⎦ −0.19315 0.21460 1.50000 0.21460 0.21460 −0.19315 0.21460 1.50000 Since this is a Dirichlet problem, A = G and, using the known values of u,we find 4 1 2 3  T b = ; −1.13192 1.13201 1.13201 −1.13192 the solution is 4 1 2 3  T Q = −0.6684 0.6684 0.6684 −0.6684 .The boundary element method 257 For the internal potential value, we have ˘ ˘ ˘ ˘ H =0.24965 = H = H = H , 11 12 13 14 ˘ ˘ ˘ ˘ G =0.08928 = G = G = G , 11 12 13 14   1 1 so that , ≈ 1.498. 2 2 The boundary values of q are significantly different from the exact values for q. The exact value for the potential is u=1+ x,which gives q(x,0) = 0,q(1,y)=1,q(x,1) = 0,q(0,y)=−1. The problem arises because, at each node, the direction of q depends on the element under consideration. This is an example of the so-called ‘corner-node problem’, and there have been a variety of ways to overcome it. We shall describe just one of these ways. The approximation for u and q as given by eqn (7.12) has the same set of basis functions for both variables; this does not have to be the case. Indeed, we may even choose different sets of nodes for the two variables. Suppose we consider the linear element and approximate u in the same manner as in Exercise 2.8. This will yield the n× n matrix H associated with the n× 1 column vector U. Now, for the q approximation, suppose we consider the set-up shown in Fig. 7.7 so that we have discontinuous elements for q. Using the notation of Exercise 7.9, the contributions to the boundary element equations in this case would be − + g q + g q . ij−1 ij j j N.B. In Exercise 7.9 we have just one value, q , and the contribution is (g + j ij−1 g )q . ij j So, now we have an n× 2n matrix G associated with a 2n× 1 column vector Q. Our overall system is again of the form HU + GQ = 0. However, care is needed in applying the boundary conditions; we now have 2n unknowns but only n conditions, either Dirichlet or Neumann. We shall not give the details here as to how we handle this; an excellent description has been given – + q q j j j j + 1 j Fig. 7.7 Nodal fluxes in adjacent elements.258 The Finite Element Method by Par´ ıs and Canas ˜ (1997), who also explain some of the other techniques for handling the corner problem. These authors also give very good descriptions of a variety of boundary elements, including isoparametric quadratics and circular elements. We have considered only Dirichlet or Neumann conditions. In principle, a Robin condition is not difficult to incorporate; see Exercise 7.10. 2 Finally, let us consider briefly Poisson’s equation ∇ u = f, which leads to  ∗ the integral equation (7.9), comprising the domain integral u fdA.If f can D 2 be written as∇ F for some function F, then Poisson’s equation can be written as Laplace’s equation 2 ∇ w=0, where w = u− F. The boundary conditions for u yield boundary conditions for w, and hence we can use the boundary element method for the solution. It is clear from this chapter that there are many similarities between the development of the system matrices for the finite element method and for the boundary element method. The major difference is that the boundary element method requires the evaluation of singular integrals. In this section, we have seen how to evaluate them for the constant element. There are other approaches, and we shall mention just two. The logarithmic integrals involved in the evaluation of G (see eqn (7.16) ii and the solution to Exercise 7.4) have been evaluated analytically. However, there is an alternative approach using a form of logarithmic Gauss quadrature given by  n g 1 f(ξ)lnξdξ≈− w f(ξ ), g g 0 g=1 where the Gauss points and associated weights are given in Appendix D. An alternative approach was developed by Telles (1987), in which a coor- dinate transformation is developed in such a way that the Jacobian vanishes at the singularity. The effect is to render the integral amenable to standard Gauss quadrature. The major advantage of this technique is that it is also applicable to the non-singular integrals, and so there is no need to distinguish between the cases in which the base node is or is not in the target element. It is particularly useful for situations where the integrands involve functions other than logarithms, for example for those integrals involving the modified Bessel function in the following section.The boundary element method 259 7.5 Time-dependent problems In this section, we give a brief introduction to the solution of the diffusion problem 1 ∂u 2 (7.28) ∇ u = in D α ∂t subject to the boundary conditions (7.29) u = g(s)on C , 1 (7.30) q = h(s)on C 2 and the initial condition (7.31) u(x,y,0) = u (x,y). 0 There are three possible approaches to dealing with such problems. Finite difference method An explicit finite difference approach in t (Smith 1985) will yield an elliptic problem of the form (k+1) (k) u − u 2 (k+1) ∇ u = α ,k =1,2,.... Δt This requires the treatment of a domain integral at each time step. We shall not take this further, but refer the reader to Wrobel (2002). Time-dependent fundamental solution It can be shown (Carslaw and Jaeger 1986) that the fundamental solution of the 2 differential operator∇ − (1/α)(∂/∂t)isgiven by   2 1 −R ∗ u (R,τ)= exp with τ = T − t, 4πατ 4ατ where the solution is sought at time T. The corresponding boundary integral equation is of the form (Kythe 1995)    t ∗ ∗ ∗ u(s,t)= α (u q− q u) dτ ds + u u dA, 0 0 C 0 D ∗ ∗ where u = u (R,T). A time-stepping scheme may now be used to develop the 0 solution (Becker 1992).260 The Finite Element Method The Laplace transform We shall define the Laplace transform as in Section 5.5 and, applying it to the problem defined by eqns (7.28)–(7.31), we obtain λ 2 ∇ u ¯− u ¯ = u in D 0 α subject to u ¯=¯ g on C , 1 ¯ q ¯ = h on C . 2 2 We shall consider the case in which u ≡ 0; if u =∇ F, we can follow the 0 0 2 approach outlined in Section 7.4. We also write p = λ/α We seek the solution of 2 2 (7.32) ∇ u ¯− p u ¯=0 in D subject to boundary conditions on C. Equation (7.32) is the modified Helmholtz equation. In Exercise 7.2, we show that the fundamental solution is 1 ∗ u ¯ = K (pR), 0 2π where K is the modified Bessel function of the second kind. 0 Since (d/dx)(K (x)) =−K (x) (Abramowitz and Stegun 1972), it follows 0 1 that 1 1 ∗ q ¯ =− pK (pR) R· n, 1 2π R and the boundary integral equation takes the form    1 u ¯ = K (pR)¯ q + pK (pR) R· nu ¯ ds. P 0 1 R C We can set up a boundary element solution, leading to a system of equations of the form ¯ ¯ HU + GQ = 0 ¯ for the boundary values of the transformed variables. Internal values U may I also be computed. Finally, the inverse transform using Stehfest’s method (see Section 5.3) will yield the solution vectors U, Q and U . IThe boundary element method 261 7.6 Exercises and solutions Exercise 7.1 The fundamental solution for the operator L is defined as that ∗ function u which satisfies ∗ Lu (r,ρ)=−δ(r,ρ) over the whole space D. The point r is fixed, the equation is referred to this position and δ(r,ρ) is the Dirac delta function, with the property  δ(r,ρ)=0, r = ρ, f(ρ) δ(r,ρ) dV (ρ)= f(r). D 2 Obtain the fundamental solution for∇ in two and three dimensions. Exercise 7.2 Obtain the fundamental solution for the modified Helmholtz equa- tion in two dimensions, 2 2 ∇ u− p u=0. Exercise 7.3 Show that the integral equation (7.7) becomes  ∗ c u = (uq − uq∗) ds, P P C where ⎧ 1, if P is inside the boundary, ⎨ c = α /2π, if P is on the boundary, P P ⎩ 0, if P is outside the boundary. Exercise 7.4 Obtain an analytical expression for the coefficient G for the ii constant element. Exercise 7.5 (Par´ ıs and Canas ˜ (1997)) In this exercise, we obtain analytic expressions for H and G (i = j) for the constant element. ij ij With reference to Fig 7.4, the angle θ varies from the value θ to θ at ij 1 2 ends 1 and 2, respectively, of the element. Obtain analytic expressions for H ij and G . ij Using these expressions, calculate the values of H and G for the problem ij ij in Example 7.1.   1 3 Exercise 7.6 Obtain the solution at , for the problem of Example 7.1. 2 4 Exercise 7.7 By considering the constant solution u≡ 1 to Laplace’s equation ˘ in a region D, show that the coefficients H for a constant element satisfy kj n ˘ H ≈ 1, kj j=1262 The Finite Element Method and use this result to check the solution to Exercise 7.6. Exercise 7.8 Use the results of Example 7.1 to solve Laplace’s equation in the unit square subject to the boundary conditions u(0,y)= u(x,0) = 0, q(1,y)=yq(x,1) = x.   1 1 Find u , . 2 2 Compare the results with the exact solution u(x,y)=1+ xy. Exercise 7.9 Obtain the coefficient matrices for a linear element. Exercise 7.10 Show how a Robin boundary condition of the form given in Section 2.3, ∂u + σ(s)u = h(s), ∂n may be incorporated into the boundary element method. Solution 7.1 We seek a solution which is dependent only on the distance from our source point. Without loss of generality, suppose that this is the origin, i.e. ∗ we seek a solution u (r) which satisfies   2 d d 2 ∗ r + ar u =−δ(r), 2 dr dr where  2 in three dimensions, a = 1 in two dimensions. First we seek a solution to the homogeneous equation which has a singularity at r =0: 2 ∗ ∗ d u a du + =0 2 dr r dr and  −1 k r in three dimensions, 3 ∗ u (r)= k lnr in two dimensions, 2 where we have ignored the additive constant in each case. Now we consider the ∗ equation for u , 2 ∗ ∇ u =−δ(r,ρ) and consider first the three-dimensional case.The boundary element method 263 Suppose we construct a small sphere V , with a surface S of radius  and   centre r; then we may write   2 ∗ ∇ u dV = −δ(r,ρ) dV V V   =−1. We now apply the divergence theorem to the integral on the left-hand side to obtain  ∗ gradu · n ˆ dS, S  so that  ∗ ∂u (r,ρ) dS =−1for r∈ V .  ∂n S  Now we can write ˆ ρ = r + RR, ˆ where R is the usual spherical polar unit vector in the radial direction relative to an origin at r. Hence, on S , R =  and ∂/∂n≡ ∂/∂R,sothat     ∗ ∗  ∂u ∂u  −1= dS = dS  ∂n ∂R S S   R= k 3 2 =− 4π 2  and 1 k = . 3 4π Consequently, in three dimensions, 1 ∗ u (r)= . 4πr In a similar manner, for the two-dimensional case we find 1 k =− 2 2π and 1 ∗ u (R)=− lnR. 2π