paidia pote perimenoume na bgei kseroume ? | Capitulo 1×04 | Shikioriori

The boundary element method for nonlinear problems

The boundary element method for nonlinear problems

Engineering Analysis with Boundary Elements 23 (1999) 365–373 The boundary element method for nonlinear problems 夽 J.T. Katsikadelis*, M.S. Nerantzak...

233KB Sizes 3 Downloads 46 Views

Engineering Analysis with Boundary Elements 23 (1999) 365–373

The boundary element method for nonlinear problems 夽 J.T. Katsikadelis*, M.S. Nerantzaki Department of Civil Engineering, National Technical University, Zografou Campus, Athens GR-157 73, Greece

Abstract In this paper a boundary-only boundary element method (BEM) is developed for solving nonlinear problems. The presented method is based on the analog equation method (AEM). According to this method the nonlinear governing equation is replaced by an equivalent nonhomogeneous linear one with known fundamental solution and under the same boundary conditions. The solution of the substitute equation is obtained as a sum of the homogeneous solution and a particular one of the nonhomogeneous. The nonhomogeneous term, which is an unknown fictitious domain source distribution, is approximated by a truncated series of radial base functions. Then, using BEM the field function and its derivatives involved in the governing equation are expressed in terms of the unknown series coefficients, which are established by collocating the equation at discrete points in the interior of the domain. Thus, the presented method becomes a boundaryonly method in the sense that only boundary discretization is required. The additional collocation points inside the domain do not spoil the pure BEM character of the method. Numerical results for certain classical nonlinear problems are presented, which validate the effectiveness and the accuracy of the proposed method. 䉷 1999 Elsevier Science Ltd. All rights reserved. Keywords: Boundary element; Nonlinear problems; Analog equation method

1. Introduction The boundary element method (BEM) has been criticized as not being capable of solving nonlinear problems. This is one of the reasons that many investigators are reluctant to be involved with BEM and use it as a computational tool, although this method has been proven to be an efficient alternative to the so called domain type methods, such as FDM and FEM, when linear problems are encountered. Effort to develop BEM methods for nonlinear problems has been given by many BEM investigators. Almost all of these methods have not avoided domain discretization. The only method that can be considered as boundary-only method is the dual reciprocity method (DRM) [1]. The term ‘‘boundary-only’’ is used in the sense that only boundary discretization is required, although collocation points inside the domain may be used to improve the solution. However, DRM works when for a non-standard linear partial differential equation or a nonlinear one it is possible to extract a standard linear partial differential operator L(·) and lump the remainder to the right-hand-side as a bodyforce term: L…u† ˆ b…x; y; u; ux ; uy ; uxx ; uxy ; uyy †

…1†

夽 This paper is dedicated to Prof. Franz Ziegler on the occasion of his 60th birthday. * Corresponding author. Tel.: ⫹30-772-1654; fax: ⫹30-772-1655. E-mail address: [email protected] (J.T. Katsikadelis)

where b(·) is, in general, a nonlinear function of its arguments. Further, DRM can be employed if the fundamental solution of the adjoint differential equation can be established, namely, a partial singular solution of the equation L*…u*† ˆ d…P ⫺ Q†

…2†

where L* (·) is the adjoint operator to L(·) and d (P ⫺ Q) is the Dirac delta function. On the basis of the aforementioned, it is apparent that the DRM cannot be employed when (a) The differential equation cannot be put in the form of the Eq. (1), e.g., uxx uyy ⫺ u2xy ˆ f …x; y† (b) The fundamental solution of Eq. (2) is not available, e.g. when the operator L* (·) has variable coefficients. Moreover, different DRM formulations and consequently different computer programs are required for different bodyforce terms as well as for different operators L*, even the order of the equation is the same. In this paper a boundary-only method is presented for solving nonlinear problems. The method is alleviated from the restrictions characterizing DRM. Simple fundamental solutions are used which depend only on the order of the

0955-7997/99/$ - see front matter 䉷 1999 Elsevier Science Ltd. All rights reserved. PII: S0955-799 7(98)00093-9

366

J.T. Katsikadelis, M.S. Nerantzaki / Engineering Analysis with Boundary Elements 23 (1999) 365–373

differential equation, e.g. for second order differential equations the fundamental solution of the Laplace equation is employed. The presented method is actually a development of the analog equation method (AEM) [2] from a domain/ boundary element method (D/BEM) to a boundary-only method (BEM). The method is illustrated by applying it to boundary value problems governed by second order nonlinear partial differential equations. However, it can also be employed to higher order equations of elliptic or hyperbolic type. Numerical results for certain problems described by strongly nonlinear differential equations were obtained to demonstrate the efficiency and accuracy of the method. Among them the minimal surface problem (Plateau’s problem) and the problem of determing a surface passing through space curve and having constant mean or constant Gaussian curvature.

The solution of Eq. (5) can be written as a sum of the  y† and a particular solution homogeneous solution u ˆ u…x; up ˆ up …x; y† of the nonhomogeneous equation: u ˆ u ⫹ up :

…7†

The particular solution is obtained from 72 up ˆ

M X

aj f j

…8†

jˆ1

which yields up ˆ

M X

aj u^ j

…9†

jˆ1

where u^j …j ˆ 1; 2; …M† are particular solutions of the equation 72 u^ j ˆ fj

j ˆ 1; 2; …M

…10†

2. The analog equation method as a boundary-only method

A particular solution of Eq. (10) can always be established when fj is specified. The homogeneous solution u is obtained from the boundary value problem.

For the sake of simplicity but without restricting its generality, the method is illustrated for problems described by second order differential equations. We consider the two-dimensional problem

72 u ˆ 0

N…u; ux ; uy ; uxx ; uyy ; uxy † ˆ g…x; y†

b 1 u ⫹ b2 q ˆ b3 ;

on G

in V

…4†

…5†

Eq. (5) indicates that the solution of Eq. (3) could be established by solving this linear equation under the same boundary condition, Eq. (4), if the fictitious source distribution b ˆ b…x; y† were known. The fictitious source can be established by working as follows. We assume bˆ

M X

aj f j

…6†

b1 u ⫹ b2 q ˆ b3 ⫺ @b1

M X

aj u^j ⫹ b2

jˆ1

M X

1

aj q^j A

on G:

jˆ1

…12† The boundary value problem Eqs. (11) and (12), is solved using the BEM. Thus, the integral representation of the solution u is given as Z ÿ   ds  u*q ⫺ uq* P 僆 {x; y} 僆 V 傼 G cu…P† ˆ⫺ G

…13† in which u* and q* are the fundamental solution of Eq. (11) and its normal derivative to the boundary, respectively, and c is the free term coefficient. Replacing u with u ⫺ up in the left hand of Eq. (13), assuming c ˆ 1 and substituting up according to Eq. (9) we obtain uˆ⫺

Z ÿ G

M X   ds ⫹ u^ j aj : u*q ⫺ uq*

…14†

jˆ1

Differentiating this equation yields ux ˆ ⫺

uy ˆ ⫺

Z ÿ G

M X   ⴱx ds ⫹ uⴱx q ⫺ uq …u^ j †x aj ;

Z 

jˆ1

where fj is a set of approximating functions and a j coefficients to be determined.

…11† 0

…3†

where N(·) is a nonlinear second order differential operator defined in a 2-D domain V describing the response of a system; bi ˆ bi …x; y†, are functions specified on boundary G; q ˆ un is the outward normal derivative of u on G. The domain V may be multiply connected. The boundary value problem, Eqs. (3) and (4), is solved using the concept of the analog equation [2], which for the problem at hand is applied as follows. Let u ˆ u…x; y† be the sought solution to the boundary value problem Eqs. (3) and (4). This function is two times continuously differentiable in V. Thus, if the Laplacian operator is applied to this function, we have 72 u ˆ b…x; y†:

in V;

uxx ˆ ⫺

G

M  X  ⴱy ds ⫹ uⴱy q ⫺ uq …u^ j †y aj ;

…15b†

jˆ1

M Z ÿ X   ⴱxx ds ⫹ uⴱxx q ⫺ uq …u^ j †xx aj ;

G

…15a†

jˆ1

jˆ1

…15c†

J.T. Katsikadelis, M.S. Nerantzaki / Engineering Analysis with Boundary Elements 23 (1999) 365–373

uyy ˆ ⫺

uxy ˆ ⫺

Z  G

Z 

M  X  ⴱyy ds ⫹ uⴱyy q ⫺ uq …u^ j †yy aj ;

…15d†

jˆ1

uⴱxy q

 ⴱxy uq



M X

367

where

  ^ : ^ ⫹ ‰b2 Š‰QŠ ‰TŠ ˆ ⫺ ‰b1 Š‰UŠ

…23†

Solving Eq. (22) we obtain …u^ j †xy aj :

…15e†

 ˆ ‰Su Š{a} ⫹ {du }; {u}

…24a†

Finally, applying Eq. (3) to M discrete points inside V, we obtain a set of M functions

 ˆ ‰Sq Š{a} ⫹ {dq } {q}

…24b†

G



ds ⫹

jˆ1

Ni …ui ; uix ; uiy ; uixx ; uiyy ; uixy † ˆ gi

i ˆ 1; 2; …M:

…16†

Using Eqs. (14) and (15a)–(15e) to evaluate u i and its derivatives at points i ˆ 1,2,…M and substituting them into Eq. (16) yields M equations of the form Ni …aj † ˆ gi

i; j ˆ 1; 2; …M

…17†

from which the coefficients a j can be established.

3. Numerical implementation of AEM The BEM with constant elements is used to approximate the boundary integrals in Eqs. (13) and (15a)–(15e). If N is the number of the boundary nodal points then Eq. (13) yields after collocating at these points ~ u}  ⫺ ‰GŠ{q}  ˆ0 ‰HŠ{

…18†

where ~ ˆ ‰HŠ ⫺ ‰IŠ{c} ‰HŠ

…19†

[H] and [G] are N × N matrices originating from the integration of the kernels on the boundary elements, {c} is a vector including the values of the coefficient c i and [I] the unit matrix. The boundary condition Eq. (12), when applied to the boundary nodal points, yields 2 3 M M X X aj u^ ij ⫹ …b2 †i aj q^ij 5 …b1 †i u i ⫹ …b2 †i q i ˆ …b3 †i ⫺ 4…b1 †i jˆ1

jˆ1

…20† or using the matrix notation



 ^ {a} ^ ⫹ ‰b2 Š‰QŠ  ⫹ ‰b2 Š{q}  ˆ {b3 } ⫺ ‰b1 Š‰UŠ ‰b1 Š{u} …21†

^ ˆ ‰q^ ij Š are N × M matrices; ^ ˆ ‰u^ij Š; ‰QŠ in which ‰UŠ ‰b1 Š; ‰b2 Š are N × N diagonal matrices and {a } the vector of the coefficients to be determined.  and Eqs. (18) and (21) may be combined to express {u}  in terms of {a }. Thus, we may write {q} " #( # ) " ( ) ~  {u} ‰0Š {0} ‰HŠ ⫺‰GŠ ˆ { a} ⫹ …22†  {q} ‰TŠ { b3 } ‰ b1 Š ‰ b2 Š

in which [Su], [Sq] are known N × M rectangular matrices and {du}, {dq} known vectors. Eqs. (14) and (15a)–(15e) when discretized and applied to M nodal points inside V give ^ a};  ⫺ ‰GŠ{q}  ⫹ ‰UŠ{ {u} ˆ ‰HŠ{u}

…25a†

 ⫺ ‰Gx Š{q}  ⫹ ‰U^ x Š{a}; {ux } ˆ ‰Hx Š{u}

…25b†

 ⫺ ‰Gy Š{q}  ⫹ ‰U^ y Š{a}; {uy } ˆ ‰Hy Š{u}

…25c†

 ⫺ ‰Gxx Š{q}  ⫹ ‰U^ xx Š{a}; {uxx } ˆ ‰Hxx Š{u}

…25d†

 ⫺ ‰Gyy Š{q}  ⫹ ‰U^ yy Š{a}; {uyy } ˆ ‰Hyy Š{u}

…25e†

 ⫺ ‰Gxy Š{q}  ⫹ ‰U^ xy Š{a}; {uxy } ˆ ‰Hxy Š{u}

…25f†

in which [G], [H], [Gx], [Hx],…,[Hxy] are known square matrices having dimensions M × N and originating from the integration of the kernel functions u* and q* and their respective derivatives on the boundary elements; ^ U^ x ],…[U^ xy ] are known matrices having dimensions [U],[ M × M, the elements of which result from the functions u^j and their derivatives. The derivatives of the kernels as well as of the functions u^ j are given in the appendix. Substituting Eqs. (24a) and (24b) into Eqs. (25a)–(25f) yields {u} ˆ ‰WŠ{a} ⫹ {w};

…26a†

{ux } ˆ ‰Wx Š{a} ⫹ {wx };

…26b†

{uy } ˆ ‰Wy Š{a} ⫹ {wy };

…26c†

{uxx } ˆ ‰Wxx Š{a} ⫹ {wxx };

…26d†

{uyy } ˆ ‰Wyy Š{a} ⫹ {wyy };

…26e†

{uxy } ˆ ‰Wxy Š{a} ⫹ {wxy };

…26f†

where ^ ‰WŠ ˆ ‰HŠ‰Su Š ⫺ ‰GŠ‰Sq Š ⫹ ‰UŠ;

…27a†

‰Wx Š ˆ ‰Hx Š‰Su Š ⫺ ‰Gx Š‰Sq Š ⫹ ‰U^ x Š;

…27b†

‰Wy Š ˆ ‰Hy Š‰Su Š ⫺ ‰Gy Š‰Sq Š ⫹ ‰U^ y Š;

…27c†

‰Wxx Š ˆ ‰Hxx Š‰Su Š ⫺ ‰Gxx Š‰Sq Š ⫹ ‰U^ xx Š;

…27d†

368

J.T. Katsikadelis, M.S. Nerantzaki / Engineering Analysis with Boundary Elements 23 (1999) 365–373

Table 1 Temperature in a square plate with temperature dependent conductivity Position x

0.1 0.3 0.5 0.7 0.9

u y

0.5 0.5 0.5 0.5 0.5

DRM

314.15 338.34 358.49 376.27 392.43

Kirchhoff

314.0 337.82 358.11 376.08 392.36

AEM M ˆ 25

M ˆ 69

M ˆ 225

314.44 339.23 359.51 376.62 392.24

314.17 337.76 358.11 376.68 392.40

314.19 338.83 359.00 375.97 391.93

While global approximating functions are independent of the collocation nodes (xj ; yj ), the RBF’s do depend on the collocation point and they appear as two-point functions. Noting that even powers of r are not RBS’s [4] we have: fj ˆ 1 ⫹ r ⫹ r3 ⫹ … ⫹ r…2m⫹1† ; rˆ



x⫺x

 j 2

…31†

  1=2 j 2 ⫹ y⫺y

‰Wyy Š ˆ ‰Hyy Š‰Su Š ⫺ ‰Gyy Š‰Sq Š ⫹ ‰U^ yy Š;

…27e†

with xj ; yj being the collocation point and x, y the point which varies during differentiation or integration. When fj are given by Eq. (31), then the functions u^ j obtained after integration are

‰Wxy Š ˆ ‰Hxy Š‰Su Š ⫺ ‰Gxy Š‰Sq Š ⫹ ‰U^ xy Š;

…27f†

u^ j ˆ

and {w} ˆ ‰HŠ{du } ⫺ ‰GŠ{dq };

…28a†

{wx } ˆ ‰Hx Š{du } ⫺ ‰Gx Š{dq };

…28b†

{wy } ˆ ‰Hy Š{du } ⫺ ‰Gy Š{dq };

…28c†

{wxx } ˆ ‰Hxx Š{du } ⫺ ‰Gxx Š{dq };

…28d†

{wyy } ˆ ‰Hyy Š{du } ⫺ ‰Gyy Š{dq };

…28e†

r2 r3 r5 r2m⫹3 ⫹ ⫹ ⫹…⫹ : 4 9 25 …2m ⫹ 3†2

…32†

4. Numerical results

{wxy } ˆ ‰Hxy Š{du } ⫺ ‰Gxy Š{dq }:

…28f†

The final step of AEM is to apply Eq. (3) to the M nodal points inside the domain. This yields the following set of M equations. o n  …29† N {u}; {ux }; {uy }; {uxx }; {uyy }; {uxy } ˆ {g}: Substituting Eqs. (26a)–(26f) into Eq. (29) yields the following M nonlinear algebraic equations including only the vector {a }, i.e., fN …{a}†g ˆ {g}:

…30†

Eq. (30) is solved numerically using any appropriate method for systems of nonlinear algebraic equations to determine the vector {a }. Once the coefficients {a } are known the solution u and its derivatives at the M nodal points are evaluated using Eqs. (26a)–(26f). For points P inside V noncoinciding with the collocation points, u and its derivatives can be evaluated from the discretized counterparts of Eqs. (14) and (15a)– (15e). The effectiveness and accuracy of the method depends on the choice of the approximating shape functions fj. Global shape functions, such as Langrange polynomials, Fourier sine and cosine series, or locally distributed functions, such as radial base functions (RBF’s) of polynomial type and thin plate spline (TPS) may be used [3]. In the development of the present method RBF’s have been utilized.

Several example problems were solved using the analytical and numerical procedure presented in the previous sections. The obtained numerical results demonstrate the efficiency, accuracy and the reliability of the presented method. In all examples the three term approximation function was employed f ˆ 1 ⫹ r ⫹ r3 :

…33†

Computationally, the method is found very efficient. All numerical results were obtained using the Fortran Power Station 4.0 of Microsoft on a Pentium 200 PC. The numerical solution of the nonlinear system of Eqs. (30) is obtained using the subroutine DNEQNF of MSIMSL. The employed algorithm is a variation of Newton’s method which uses a finite difference approximation to the Jacobian. 4.1. Example 1: Heat flow in bodies with nonlinear material properties In practical heat transfer problems the conductivity may be a function of the temperature, k ˆ k…u†. In this case the steady state heat transfer equation, neglecting the heat sources, can be written as ÿ    …34† kux x ⫹ kuy ˆ 0: y

Assuming that k varies according to the law   k ˆ ko 1 ⫹ b…u ⫺ uo †=uo

…35†

where uo is the initial temperature and ko, b are material constants, Eq. (34) becomes  k  …36† k72 u ⫹ o b u2x ⫹ u2y ˆ 0: uo The nonlinear Eq. (36) is solved for a square plane body with unit side length 0 ⱕ x, y ⱕ 1 under mixed boundary

J.T. Katsikadelis, M.S. Nerantzaki / Engineering Analysis with Boundary Elements 23 (1999) 365–373 Table 2 Numerical results of the field function u and its derivatives in Example 2. Upper row: Computed: Lower row: Exact x

u

ux

uy

0.2273

6.610 6.610 6.579 6.579 6.516 6.516 6.421 6.420 6.292 6.290 6.127 6.124 5.917 5.917 5.655 5.668 5.339 5.369 4.979 5.010 4.568 4.579

⫺ 0.034 ⫺ 0.034 ⫺ 0102 ⫺ 0.103 ⫺ 0.174 ⫺ 0.174 ⫺ 0.246 ⫺ 0.247 ⫺ 0.321 ⫺ 0.325 ⫺ 0.408 ⫺ 0.408 ⫺ 0.518 ⫺ 0.499 ⫺ 0.638 ⫺ 0.601 ⫺ 0.745 ⫺ 0.719 ⫺ 0.838 ⫺ 0.861 ⫺ 0.987 ⫺ 1.042

⫺ 0.367 ⫺ 0.378 ⫺ 0.381 ⫺ 0.380 ⫺ 0.382 ⫺ 0.383 ⫺ 0.385 ⫺ 0.389 ⫺ 0.392 ⫺ 0.397 ⫺ 0.408 ⫺ 0.408 ⫺ 0.436 ⫺ 0.422 ⫺ 0.469 ⫺ 0.441 ⫺ 0.500 ⫺ 0.465 ⫺ 0.523 ⫺ 0.499 ⫺ 0.554 ⫺ 0.545

0.6818 1.1364 1.5909 2.0455 2.5000 2.9545 3.4091 3.8636 4.3182 4.7727

uxx ⫺ 0.177 ⫺ 0.151 ⫺ 0.147 ⫺ 0.153 ⫺ 0.158 ⫺ 0.155 ⫺ 0.158 ⫺ 0.165 ⫺ 0.168 ⫺ 0.175 ⫺ 0.193 ⫺ 0.190 ⫺ 0.224 ⫺ 0.211 ⫺ 0.239 ⫺ 0.240 ⫺ 0.230 ⫺ 0.282 ⫺ 0.262 ⫺ 0.347 ⫺ 0.465 ⫺ 0.455

uyy ⫺ 0.151 ⫺ 0.172 ⫺ 0.180 ⫺ 0.173 ⫺ 0.179 ⫺ 0.176 ⫺ 0.185 ⫺ 0.179 ⫺ 0.190 ⫺ 0.184 ⫺ 0.193 ⫺ 0.190 ⫺ 0.206 ⫺ 0.199 ⫺ 0.241 ⫺ 0.210 ⫺ 0.283 ⫺ 0.226 ⫺ 0.290 ⫺ 0.249 ⫺ 0.269 ⫺ 0.283

conditions u…0; y† ˆ 300;

un …x; 0† ˆ 0;

u…1; y† ˆ 400;

un …x; 1† ˆ 0

…37†

and for k0 ˆ 1, b ˆ 3, u0 ˆ 300. The numerical results are presented in Table 1 as compared with those obtained using the DRM and Kirchhoff’s transformation [1]. The results were obtained using 100 constant boundary elements and M ˆ 25, 69 and 225 internal collocation points.

4.2. Example 2: Determination of a surface with constant mean curvature The determination of a surface u ˆ u…x; y† bounded by one or more nonintersecting space curves and having constant mean curvature k is obtained from the following

369

boundary value problem [5]    3=2 …1 ⫹ u2y †uxx ⫺ 2ux uy uxy ⫹ 1 ⫹ u2x uyy ⫺ k 1 ⫹ u2x ⫹ u2y ˆ0

in V;

…38a†

u ˆ a3

on G:

…38b†

This equation is of elliptic type, because its discriminant, namely,    1 ⫹ u2x 1 ⫹ u2y ⫺ u2x u2y ˆ 1 ⫹ u2x ⫹ u2y is greater than zero. Eq. (38a) is solved when V is the square 0 ⱕ x, y ⱕ 5 under the boundary conditions:  1=2  1=2 u…5; y† ˆ 25 ⫺ y2 ; u…0; y† ˆ 50 ⫺ y2 ; …39†     2 1=2 2 1=2 ; u…x; 5† ˆ 25 ⫺ x u…x; 0† ˆ 50 ⫺ x p and mean curvature k ˆ ⫺ 2=5. The exact solution is  i1=2 h …40† u ˆ 50 ⫺ x2 ⫹ y2 which represents a spherepwith  center the origin of the coordinates and radius R ˆ 5 2: Numerical results for u and its derivatives on the central line y ˆ 2.5 are presented in Table 2. They have been obtained using 100 constant boundary elements and 121 interior points uniformly distributed in the square domain (xi ˆ Dx=2 ⫹ …i ⫺ 1†Dx; yj ˆ Dy=2 ⫹ …j ⫺ i†Dy; i ˆ 1; 2; …; 11; Dx ˆ Dy ˆ 5=11). The accuracy for u is very good (the maximum error does not exceeds 0.5%). The accuracy of the computed derivatives decreases as their order increases. Moreover, in Table 3 the dependence of u and its derivatives computed at the center of the domain (x ˆ y ˆ 2.5) on the number of the interior collocation points are presented. The principal curvature is evaluated from the expression uxx k1 ˆ k2 ˆ kx ˆ ky ˆ ÿ …41† 1=2 :  2 1 ⫹ ux 1 ⫹ u2x ⫹ u2y The results were obtained using 100 boundary elements. Increase of the number of the boundary elements did not improve the accuracy of the results.

Table 3 Dependence of u and its derivatives at x ˆ y ˆ 2.5 on the number M of the interior collocation points M 25 u ux ˆ uy uxx ˆ uyy uxy k1 ˆ k2

6.1096 ⫺ 0.42502 ⫺ 0.19416 ⫺ 0.02556 ⫺ 0.14095

49 6.1182 ⫺ 0.42660 ⫺ 0.19736 ⫺ 0.04393 ⫺ 0.14296

81

121

169

Exact

6.1248 ⫺ 0.41869 ⫺ 0.19645 ⫺ 0.50840 ⫺ 0.14382

6.1276 ⫺ 0.40879 ⫺ 0.19361 ⫺ 0.04795 ⫺ 0.14361

6.1272 ⫺ 0.39884 ⫺ 0.18890 ⫺ 0.03099 ⫺ 0.14195

6.1237 ⫺ 0.40825 ⫺ 0.19052 ⫺ 0.02721 ⫺ 0.14142

370

J.T. Katsikadelis, M.S. Nerantzaki / Engineering Analysis with Boundary Elements 23 (1999) 365–373

Table 4 Numerical results for minimal surface and its derivatives along the line y ˆ 5-x x

u

ux

uy

uxx

uyy

uxy

4.7727 4.3182 3.8636 3.4091 2.9545 2.50 2.0455 1.5909 1.1364 0.6818 0.2273

22.492 16.825 11.195 6.4336 2.8133 0.3 × 10 ⫺13 ⫺ 2.8133 ⫺ 6.4336 ⫺ 11.195 ⫺ 16.825 ⫺ 22.492

10.334 9.7070 8.0194 5.9094 4.0202 2.9247 2.8239 3.3019 3.6216 3.0993 1.4114

⫺ 1.4114 ⫺ 3.0933 ⫺ 3.6217 ⫺ 3.3019 ⫺ 2.8239 ⫺ 2.9247 ⫺ 4.0202 ⫺ 5.9094 ⫺ 8.0194 ⫺ 9.7070 ⫺ 10.334

1.1235 4.0170 3.9059 3.5538 2.6680 1.2765 0.12705 ⫺ 0.75746 0.49892 1.2752 1.2466

⫺ 1.2466 ⫺ 1.2752 ⫺ 0.49892 0.75746 ⫺ 0.12705 ⫺ 1.2765 ⫺ 2.6680 ⫺ 3.5538 ⫺ 3.9059 ⫺ 4.0170 ⫺ 1.1235

4.4914 1.3101 ⫺ 0.38827 ⫺ 1.1536 ⫺ 0.95847 0.1 × 10 ⫺12 0.95847 1.1536 0.38827 ⫺ 1.3101 ⫺ 4.4914

4.3. Example 3: The problem of minimal surface This problem concerns with the determination of a surface bounded by one or more nonintersecting space curves and having minimum area. The governing differential equation resulting as the Euler–Langrange condition for minimization of the integral  Z  1 ⫹ u2x ⫹ u2y dx dy …42† Aˆ

that of soap bubble constrained between bounding contours (Plateau’s problem [6]). Eq. (43a) is solved when the domain V is the square 0 ⱕ x, y ⱕ 5 under the boundary conditions: u…0; y† ˆ ⫺y2 ;

u…5; y† ˆ 25 ⫺ y2 ;

u…x; y† ˆ x2 ⫺ 25:

u…x; 0† ˆ ⫺y2 ; …44†

V

is     1 ⫹ u2y uxx ⫺ 2ux uy uxy ⫹ 1 ⫹ u2x uyy ˆ 0

in V …43a†

and u ˆ a3

on G:

…43b†

Comparing Eq. (43a) with Eq. (38a) it follows that the minimal surface can be defined as one for which the mean curvature is zero. The physical analogue of the minimal surface is

The boundary conditions are antisymmetric with respect to the diagonal x ˆ y. The numerical results were obtained using 100 constant boundary elements and 121 interior points uniformly distributed in the square domain. They are presented in graphical form in Figs. 1 and 2. One can readily see that the solution is antisymmetric with respect to the diagonal x ˆ y. Numerical results along the line y ˆ 5 ⫺ x are also presented in Table 4. Finally, the convergence of u and its derivatives at the center of the domain with increasing

Fig. 1. The surface of the soap bubble in Example 3 (Section 4.3).

J.T. Katsikadelis, M.S. Nerantzaki / Engineering Analysis with Boundary Elements 23 (1999) 365–373

371

Fig. 2. Contours of the soap bubble in Example 3 (Section 4.3).

number of the interior collocation points is shown in Table 5. 4.4. Example 4: Determination of a surface with constant Gaussian curvature The Gaussian or total curvature of a surface u ˆ u…x; y† defined as the product of the two principal curvatures is given as [7] Kˆ

f 2 ⫺ eg F 2 ⫺ EG

…45†

Substituting Eq. (46) into Eq. (45) gives Kˆ 

uxx uyy ⫺ u2xy 2 : 1 ⫹ u2x ⫹ u2y

…47†

Thus, the problem of determining a surface u(x,y) bounded by one or more nonintersecting space curves with a given constant Gaussian curvature is obtained from the following boundary value problem  2 uxx uyy ⫺ u2xy ⫺ K 1 ⫹ u2x ⫹ u2y ˆ 0

in V;

…48a†

where, E, F, G and e, f, g are the coefficients of the first and second fundamental forms, respectively, which can be expressed in terms of u as

u ˆ a3

E ˆ 1 ⫹ u2x ;

Eq. (48a) is solved when the domain V is the square 0 ⱕ x,

F ˆ ux u y ;

G ˆ 1 ⫹ u2y ;

 1=2 e ˆ uxx = 1 ⫹ u2x ⫹ u2y ;

 1=2 f ˆ uxy = 1 ⫹ u2x ⫹ u2y ;

 1=2 g ˆ uyy = 1 ⫹ u2x ⫹ u2y :

…46†

Table 5 Convergence of u and its derivatives at x ˆ y ˆ 2:5 M ˆ 25 u ux uy uxx uyy uxy

M ˆ 49 ⫺13

0.1 × 10 2.9890 ⫺ 2.9890 1.1541 1.1541 0.2 × 10 ⫺13

M ˆ 81 ⫺12

0.3 × 10 2.9738 ⫺ 2.9738 1.1914 ⫺ 1.1914 0.4 × 10 ⫺12

M ˆ 121 ⫺14

0.2 × 10 2.9488 ⫺ 2.9488 1.2646 ⫺ 1.2646 0.7 × 10 ⫺13

0.5 × 10 ⫺13 2.9247 2.9247 1.2765 ⫺ 1.2765 0.3 × 10 ⫺12

on G:

…48b†

Table 6 Numerical results for u in Example 4 (Section 4.4) x

0.25 0.75 1.25 1.75 2.25 2.75 3.25 3.75 4.25 4.75

y ˆ 1.25

y ˆ 3.75

AEM

Exact

AEM

Exact

6.92 6.82 6.70 6.56 6.40 6.21 6.01 5.71 5.39 5.03

6.95 6.91 6.84 6.73 6.58 6.39 6.15 5.86 5.51 5.08

5.94 5.81 5.68 5.55 5.39 5.18 4.91 4.59 4.20 3.68

5.98 5.94 5.86 5.73 5.55 5.32 5.03 4.67 4.22 3.65

372

J.T. Katsikadelis, M.S. Nerantzaki / Engineering Analysis with Boundary Elements 23 (1999) 365–373

h i1=2 r ˆ …x ⫺ j†2 ⫹ …y ⫺ h†2

y ⱕ 5 under the boundary conditions u…0; y† ˆ …50 ⫺ y † ;

u…5; y† ˆ …25 ⫺ y † ;

u…x; 0† ˆ …50 ⫺ x † ;

u…x; 5† ˆ …25 ⫺ x †

2 1=2 2 1=2

x; y 僆 V; j; h 僆 G …A:1†

2 1=2

and Gaussian curvature K ˆ 1/50. The exact solution is  i1=2 h u ˆ 50 ⫺ x2 ⫹ y2

2 1=2

…49†

rx ˆ

…50†

5. Conclusions A general boundary-only BEM has been developed to solve nonlinear problems in engineering. Its formulation and solution procedure is demonstrated for problems described by second order nonlinear differential equations. The method is based on the concept of the analog equation, according to which the nonlinear equation is substituted by a linear one which is solved using BEM. The main features of the method are: 1. The method is pure BEM in the sense that only boundary discretization is required. 2. The method is straight forward. That is, the formulation of the solution procedure depends only on the order of the operator and not on the specific problem described by it. 3. There is no need to extract a standard partial operator as it happens with DRM. 4. Linear problems for inhomogeneous bodies leading to equations with variable coefficients can also be treated as special cases using the same formulation. 5. In all problems the employed fundamental solution is that of the standard differential operator, e.g. for equations of second order that of Laplace. 6. Accurate results are obtained using a relatively small number of interior collocation points. The accuracy of the results depends highly on the method for solving the system of nonlinear equations. The results are improved when high precision arithmetic is used. 7. The computer program is the same for all problems.

A.1. Derivatives of the kernel functions The derivatives of the distance

x⫺j ; r

ry ˆ

y⫺n ; r

rj ˆ ⫺rx ;

rh ˆ ⫺ry ;

ry2 rx ry r2 ; ryy ˆ x ; rxy ˆ ⫺ ; r r r     rn ˆ ⫺ rx nx ⫹ ry ny ; rt ˆ ⫺ ⫺rx ny ⫹ ry nx :

rxx ˆ

which represents a sphere withpcenter the origin of the coordinates and radius R ˆ 5 2. Numerical results are presented in Table 6 as compared with the exact ones. They have been obtained using 100 constant boundary elements and 100 interior collocation points uniformly distributed in the square domain …xi ˆ Dx=2 ⫹ …i ⫺ 1†Dx; yj ˆ Dy=2 ⫹ …j ⫺ 1†Dy; i; j ˆ 1; 2; …10; Dx ˆ Dy ˆ 5=10†:

Appendix A

can be evaluated from the following relations

…A:2†

nx, ny are directional cosines of the outward normal to the boundary. Using Eqs. (A.2) the derivatives of the fundamental solution, Eq. (12), may be expressed as uⴱx ˆ

1 rx ; 2p r

uⴱy ˆ

1 ry ; 2p r

2 2 1 ry ⫺ rx ; uⴱyy ˆ ⫺uⴱxx ; 2 2p r   2 2 1 ry ⫺ rx rn ⫹ 2rx ry rt ˆ⫺ ; p r3

uⴱxx ˆ uⴱnxx

uⴱnyy ˆ ⫺uⴱnxx ;   ry2 ⫺ rx2 rt ⫺ 2rx ry rn 1 : uⴱnxy ˆ ⫺ p r3

uⴱxy ˆ ⫺

1 rx ry ; p r2

…A:3†

A.2. Derivatives of the functions u^j For the three term approximation function Eq. (32) becomes u^ j ˆ

r2 r3 r5 ⫹ ⫹ 4 9 25

…A:4†

which gives by differentiation    1 1 1 ⫹ r ⫹ r 2 x ⫺ xj ; u^ jx ˆ 2 3 4    1 1 1 ⫹ r ⫹ r 2 y ⫺ yj ; u^ jy ˆ 2 3 4   1 1 2 1 3 4 j ⫹ ⫹ r ⫹ u^ xx ˆ 2 r2 3 r 4 5   2  1 1  1 1 1 1  j j 2 ⫹ ⫹ r y ⫺ y ⫹ ;  x⫺x ⫹ 2 r2 3 r 4 5  u^ jyy ˆ

 1 1 1 1 1 1 ⫹ ⫹ r ⫹ 2 r2 3 r 4 5   2  1 1  2 1 3 4  j 2 ⫹ ⫹ r y ⫺ y ⫹ ;  x ⫺ xj ⫹ 2 r2 3 r 4 5

J.T. Katsikadelis, M.S. Nerantzaki / Engineering Analysis with Boundary Elements 23 (1999) 365–373

 u^ jxy ˆ

1 1 1 3 ⫹ ⫹ 3 r 2 5



  x ⫺ xj y ⫺ yj :

…A:5† [2]

It can be readily proved that lim u^ jx ˆ 0; r!0

lim u^ jyy ˆ r!0

1 ; 2

lim u^ jy ˆ 0; r!0

lim u^ jxx ˆ r!0

1 ; 2

lim u^jxy ˆ 0:

[3]

[4]

r!0

[5]

References [1] Partridge PW, Brebbia CA, Wrobel LC. The dual reciprocity boundary

[6] [7]

373

element method. Southampton: Computational Mechanics Publications, 1992. Katsikadelis JT. The analog equation method – A powerful BEM-based solution technique for solving linear and nonlinear engineering problems. In: Brebbia CA, editor. Boundary Element Method XVI, Southampton: CLM Publications, 1994. pp. 167. Partridge PW. Approximation functions in the dual reciprocity method. International Journal of Boundary Elements Communications 1997;8:1–4. Yamada T, Wrobel LC, Power H. On the convergence of the dual reciprocity method. Engineering Analysis with Boundary Elements 1994;13:291–298. Davis HT. Introduction to nonlinear differential and integral equations. New York: Dover Publications Inc, 1962 Chap. 14. Rado T. On the problem of plateau. Berlin: Springer-Verlag, 1932. Lass H. Vector and tensor analysis. New York: McGraw-Hill, 1950 Chap. 3.