International Scholarly Research Notices

International Scholarly Research Notices / 2011 / Article

Research Article | Open Access

Volume 2011 |Article ID 295904 |

V. V. Zozulya, "Numerical Solution of the Kirchhoff Plate Bending Problem with BEM", International Scholarly Research Notices, vol. 2011, Article ID 295904, 14 pages, 2011.

Numerical Solution of the Kirchhoff Plate Bending Problem with BEM

Academic Editor: T. K. Basak
Received27 Jan 2011
Accepted24 Mar 2011
Published15 Jul 2011


Direct approach based on Betty's reciprocal theorem is employed to obtain a general formulation of Kirchhoff plate bending problems in terms of the boundary integral equation (BIE) method. For spatial discretization a collocation method with linear boundary elements (BEs) is adopted. Analytical formulas for regular and divergent integrals calculation are presented. Numerical calculations that illustrate effectiveness of the proposed approach have been done.

1. Introduction

The BIE method and its numerical realization BEM are powerful tools for analysis of the wide range problems in mechanics and engineering [1, 2]. The main advantage of the BEM in comparison with other numerical methods (e.g., finite element method (FEM)) is that in the BEM only boundary of the domain and unknown functions on it have to be approximated. It leads to reduction of the problem dimension, that is, a 3D problem is transformed to a 2D one, and so forth. As result the amount of calculation is significantly reduced.

During the last decade, the BEM has been established as a robust numerical method for solution of the elastic thin plate problems [3]. Various aspects of the BEM application to the problems of thin Kirchhoff’s plates bending have been reported in (see [411, 22], etc). In these and other publications calculations the corresponding integrals have been doing with application of numerical Gaussian’s quadratures. It reduces accuracy and effectiveness of the BEM. In [10, 12] analytical integration has been applied, and exact formulas for the integrals evaluation have been obtained.

Another problem which gives many troubles in application of the BEM to the Kirchhoff’s plates bending is the divergent integrals with singularities of the following type: weakly singular 𝑂(ln𝑟), strongly singular 𝑂(𝑟1), and hypersingular 𝑂(𝑟2). For their calculation special methods have to be applied. Various methods for divergent integrals regularization have been developed in (see [1316], etc). Application of the regularization algorithms for the divergent integrals evaluation in the BEM solution of the plate bending problems can be found in [5, 6, 8, 17]. Information related to various aspects of the numerical solution of integral equations can be found in [14, 18, 19] and related to integral calculation numerical using various quadrature formulas in [20, 21].

In this paper, we present a general direct BEM approach for thin Kirchhoff’s plate bending problems. For evaluation of the regular integrals we use analytical formulas which give accurate results and need less time for calculation. For divergent integrals regularization we use special method developed in [15, 16]. This method is based on theory of distribution and gives simple analytical equations for evaluation of weakly singular, strongly singular, and hypersingular integrals in the same way. Numerical examples that illustrate effectiveness of the developed BEM for calculation of the thin plate bending problems are presented.

2. The Plate Equations

Let us consider thin elastic plate with thickness 2. The plate is subjected by load of the density 𝑝(𝐱), transversal to its middle surface Ω, as it is shown on Figure 1. We suppose that deflection of the plate is small in comparison with its thickness, and Kirchhoff’s hypothesis takes place [8].

We consider element that cuts out from the plate by two pairs of planes which are parallel to the planes 𝑥1𝑥3 and 𝑥2𝑥3 (Figure 2). Then equilibrium of the plate element can be expressed by the equations𝜕𝑖𝑞𝑖(𝐱)+𝜕𝑗𝑞𝑗𝜕(𝐱)+𝑝(𝐱)=0,𝑖𝑚𝑖𝑗(𝐱)𝜕𝑗𝑚𝑗(𝐱)+𝑞𝑗𝜕(𝐱)=0,𝑗𝑚𝑗𝑖(𝐱)+𝜕𝑖𝑚𝑖(𝐱)𝑞𝑖(𝐱)=0.(1) Here and further on 𝑤 is a deflection of the plate, 𝑚𝑖𝑗 are the bending and twisting moments, and 𝑞𝑖(𝐱) is a shear force. Indices 𝑖,𝑗 possess the value 1,2 and correspond to coordinates 𝑥1,𝑥2. Arrows on Figure 2 show positive directions of the bending and twisting moments.

From (1) follows inhomogeneous biharmonic differential equation for the deflection of the plate in the form [8] 𝑝ΔΔ𝑤=𝐷,𝐱Ω,(2) where ΔΔ is the biharmonic differential operator, Δ=𝜕𝑖𝜕𝑖 is the Laplace operator, 𝜕i=𝜕/𝜕𝑥𝑖 is the partial derivative with respect to coordinate 𝑥𝑖, 𝜈 is Poisson’s ratio, D=2𝐸2/3(1𝜈2) is the plate stiffness, and 𝐸 is the elastic modulus.

The plate deflection is related to other parameters of the plate by relations 𝜃𝑖=𝜕𝑖𝑤,𝑚𝑖,𝑗=𝐷(1𝜈)𝜕𝑖𝜕𝑗𝑤+𝜈𝛿𝑖𝑗𝜕𝑘𝜕𝑘𝑤,𝑞𝑖=𝐷𝜕𝑖Δ𝑤,(3) where 𝜃𝑖 is the slope of the element perpendicular to the middle surface of the plate.

Stress-strain state of the plate is defined by stress and strain tensors. Due to Kirchhoff’s hypothesis they are related to the deflection by the equations𝜎𝑖𝑗=𝑥3𝜆𝛿𝑖𝑗Δ𝑤+2𝜇𝜕𝑖𝜕𝑗𝑤,𝜀𝑖𝑗=𝑥3𝜕𝑖𝜕𝑗𝑤.(4)

Moments and shear force are expressed through stress tensor in the form 𝑚𝑖𝑗=𝜎𝑖𝑗𝑥3𝑑𝑥3,𝑞𝑖=𝜎𝑖3𝑥3𝑑𝑥3,(𝑖=1,2).(5)

Generalized shear force may be introduced by the relation 𝑣𝑛=𝑞𝑛𝜕𝑡𝑚𝑡.(6)

Correct statement of the plate bending problem supposes that (2) is supplemented by suitable conditions on the boundary 𝜕Ω. In the Kirchhoff’s theory of plates boundary conditions are usually prescribed to the functions [8]𝑤,𝜃𝑛=𝜃𝑖𝑛𝑖=𝜕𝑛𝑚𝑤,𝑛=𝑚𝑖𝑗𝑛𝑖𝑛𝑗𝜕=𝐷𝑛𝜕𝑛𝑤+𝜈𝜕𝑡𝜕𝑡𝑤,𝑣𝑛=𝑞𝑖𝑛𝑖𝜕𝑡𝑚𝑖𝑗𝑛𝑖𝑡𝑗=𝐷𝜕𝑛Δ𝑤+(1𝜈)𝜕2𝑡𝑤,(7) where 𝑛𝑖и𝑡𝑖 are the components of vectors normal and tangential to the boundary 𝜕Ω, respectively, 𝜕𝑛=𝑛𝑖𝜕𝑖 and 𝜕𝑡=𝑡𝑖𝜕𝑖 are normal and tangential derivatives, respectively.

For arbitrary part of the boundary we have (see Figure 3) 𝑚𝑖=𝑚𝑖𝑗𝑛𝑗,𝑚𝑛=𝑚𝑖𝑛𝑖=𝑚𝑖𝑗𝑛𝑖𝑛𝑗,𝑞𝑖=𝑞𝑖𝑛𝑖,𝑚𝑡=𝑚𝑖𝑡𝑖=𝑚𝑖𝑗𝑛𝑖𝑡𝑗,𝑚𝑖=𝑚𝑛𝑛𝑖+𝑚𝑡𝑡𝑖,𝜃𝑡𝑖=𝜃𝑖𝑡𝑖,𝜃𝑖=𝜃𝑛𝑛𝑖+𝜃𝑡𝑡𝑖.(8)

In the case if boundary contains corner points, corner forces 𝑞𝑐 which are generated by 𝑚𝑡 in the points where tangential vector has break have to be taken into account in the expression for 𝑣𝑛. Really from the Figure 4 follows that the 𝑚𝑡 acts on different sides of any angle 𝛼 in opposite directions and 𝑚𝑡 and 𝑚+𝑡 related to the different sides of the angle therefore, 𝑞𝑐=𝑚𝑡𝑚+𝑡.

Under the Kirchhoff assumptions boundary conditions can be defined as follows:𝑤=0,𝜃𝑛𝑞=0,clampedcontour𝑐=0,𝑤𝑐=0(9)𝑤=0,𝑚𝑛𝑞=0,simply-supportedcontour𝑐0,𝑤𝑐𝑣=0(10)𝑛=0,𝑚𝑛𝑞=0,freecontour𝑐=0,𝑤𝑐0(11)

Each of the boundary conditions (9)–(11) together with differential equation for the deflection of the plate (2) represents corresponding boundary value problem.

3. Betty’s Reciprocal Theorem and Elastic Potentials

The boundary-value problems (2), (9)–(11) can be solved by the BIE method. In order to transform a boundary-value problem to the BIE, Betty’s reciprocal theorem can be used in the following form [1, 3]:Ω𝑝𝑤𝑑Ω+𝜕Ω𝑞𝑛𝑤+𝑚𝑖𝜃𝑖=𝑑𝑆Ω𝑝𝑤𝑑Ω+𝜕Ω𝑞𝑛𝑤+𝑚𝑖𝜃𝑖𝑑𝑆.(12) In expression (12) functions with accent 𝑝,𝑤,𝜃𝑛,𝑚𝑛,𝑣𝑛 and without accent 𝑝,𝑤,𝜃𝑛,𝑚𝑛,𝑣𝑛 are parameters of two different states of the elastic plate, respectively.

Let us transform (12) to the form suitable for the case of presence mentioned above corner points. For that purpose we consider series of integral transformationsΩ𝑚𝑖𝜃𝑖𝑑𝑆=𝜕Ω𝑚𝑛𝑛𝑖+𝑚𝑡𝑡𝑖𝜃𝑛𝑛𝑖+𝜃𝑡𝑡𝑖=𝑑𝑆𝜕Ω𝑚𝑛𝜃𝑛+𝑚𝑡𝜃𝑡=𝑑𝑆𝜕Ω𝑚𝑛𝜃𝑛+𝑚𝑡𝜕𝑡𝑤𝑑𝑆.(13) Last integral in (13) can be represented in the form𝜕Ω𝑚𝑡𝜕𝑡𝑤𝑑𝑆=𝜕Ω𝜕𝑡𝑚𝑡𝑤𝑑𝑆𝜕Ω𝜕𝑡𝑚𝑡𝑤𝑑𝑆.(14) For close contour the first right integral in (14) is equal to 𝐼=𝜕Ω𝜕𝑡𝑚𝑡𝑤𝑑𝑆=𝜕Ω𝑑𝑚𝑡𝑤𝑑𝑆=0.(15) For contour with corner points 𝐼=𝑚𝑡𝑤𝑚+𝑡𝑚𝑤=𝑡𝑚+𝑡𝑤=𝑞𝑐𝑤𝑐.(16) Finally integral (13) can be presented in the form 𝜕Ω𝑚𝑖𝜃𝑖𝑑𝑆=𝜕Ω𝑚𝑛𝜃𝑛𝜕𝑡𝑚𝑡𝑤𝑑𝑆+𝑞𝑐𝑤𝑐.(17)

In the same way can be transformed integral𝜕Ω𝑚𝑖𝜃𝑖𝑑𝑆=𝜕Ω𝑚𝑛𝜃𝑛𝜕𝑡𝑚𝑡𝑤𝑑𝑆+𝑞𝑐𝑤𝑐.(18)

With taking into account representations (17) and (18) Betty’s reciprocal relation (12) can be extended for the case of the corner points presence and can be represented in the form Ω𝑝𝑤𝑑Ω+𝜕Ω𝑣𝑛𝑤+𝑚𝑛𝜃𝑛𝑑𝑆+𝑞𝑐𝑤𝑐=Ω𝑝𝑤𝑑Ω+𝜕Ω𝑣𝑛𝑤+𝑚𝑛𝜃𝑛𝑑𝑆+𝑞𝑐𝑤𝑐.(19)

From Betty’s reciprocal relation (13) one can obtain Somigliano’s type identity for the deflection of the plate 𝑤. To do that we choose as the first state of plate (without accent) that state, which has to be calculated. The second state of the plate (with accent) is auxiliary and correspond to infinite plate subjected to action of the concentrated force of unit density that applied at point 𝐱 in the direction normal to the middle surface. Functions that correspond to the second state are fundamental solutions for the plate. For convenience the following notations will be used: 𝑝=𝛿(𝐲,𝐱),𝑤=𝑊(𝐲,𝐱),𝜃𝑛𝑚=Θ(𝐲,𝐱),𝑛=𝑀(𝐲,𝐱),𝑣𝑛=𝑉(𝐲,𝐱).(20)

Taking into account that Dirac’s delta function 𝛿(𝐱,𝐲) transform regular functions in the following way: Ω𝑤(𝐲)𝛿(𝐲,𝐱)𝑑Ω=𝑤(𝐱),(21) Betty’s reciprocal relation (19) can be transformed to the Somigliano’s type identity for the plate deflection 𝑤(𝐱)=𝜕Ω𝑣𝑛(𝐲)𝑊(𝐲,𝐱)+𝑚𝑛(𝐲)Θ(𝐲,𝐱)𝜃𝑛+(𝐲)𝑀(𝐲,𝐱)𝑤(𝐲)𝑉(𝐲,𝐱)𝑑𝑆Ω𝑝+(𝐲)𝑊(𝐲,𝐱)𝑑Ω𝐾𝑘=1𝑞𝑐𝑊𝑐(𝐲,𝐱)𝑀𝑐(𝐲,𝐱)𝑤𝑐𝑘,𝐱R2𝜕Ω.(22) Here 𝐾 is the number of corner points on the contour of the plate 𝜕Ω.

The slope of the element perpendicular to the middle surface of the plate 𝜃𝑛(𝐱) can be obtained applying operator of the normal derivative 𝜕𝑛 to (22). We assume that point 𝐱𝑅2𝜕Ω and therefore kernels in (22) are sufficiently smooth and differentiation with respect to parameter 𝐱 is correct. As result we obtain 𝜃𝑛(𝐱)=𝜕Ω𝑣𝑛(𝐲)𝑊𝑛(𝐲,𝐱)+𝑚𝑛(𝐲)Θ𝑛(𝐲,𝐱)𝜃𝑛(𝐲)𝑀𝑛(𝐲,𝐱)𝑤(𝐲)𝑉𝑛+(𝐲,𝐱)𝑑𝑆Ω𝑝(𝐲)𝑊𝑛+(𝐲,𝐱)𝑑Ω𝐾𝑘=1𝑞𝑐𝑊𝑐𝑛(𝐲,𝐱)𝑀𝑐𝑛(𝐲,𝐱)𝜕𝑛𝑤𝑐𝑘,𝐱𝑅2𝜕Ω.(23) The kernels 𝑊𝑛, Θ𝑛, 𝑀𝑛, and 𝑉𝑛 can be obtained by applying operator of the normal derivative 𝜕𝑛 with respect to 𝐱 to corresponding fundamental solutions in (22).

Further in order to construct BIE limit transition to the boundary 𝐱𝜕Ω will be done.

For convenience and compactness of the BIE consideration we introduce the following contour and area potentials: 𝑊𝑣𝑛=,𝐱,𝜕Ω𝜕Ω𝑣𝑛(Θ𝑚𝐲)𝑊(𝐲,𝐱)𝑑𝑆,𝑛=,𝐱,𝜕Ω𝜕Ω𝑚𝑛(𝑀𝜃𝐲)Θ(𝐲,𝐱)𝑑𝑆,𝑛=,𝐱,𝜕Ω𝜕Ω𝜃𝑛(𝐲)𝑀(𝐲,𝐱)𝑑𝑆,𝑉(𝑤,𝐱,𝜕Ω)=𝜕Ω𝑊𝑤(𝐲)𝑉(𝐲,𝐱)𝑑𝑆,𝑛𝑣𝑛=,𝐱,𝜕Ω𝜕Ω𝑣𝑛(𝐲)𝑊𝑛Θ(𝐲,𝐱)𝑑𝑆,𝑛𝑚𝑛=,𝐱,𝜕Ω𝜕Ω𝑚𝑛(𝐲)Θ𝑛(𝐲,𝐱)𝑑𝑆,(24)𝑀𝑛𝜃𝑛=,𝐱,𝜕Ω𝜕Ω𝜃𝑛(𝐲)𝑀𝑛(𝑉𝐲,𝐱)𝑑𝑆,𝑛(𝑤,𝐱,𝜕Ω)=𝜕Ω𝑤(𝐲)𝑉𝑛(𝐲,𝐱)𝑑𝑆,𝑊(𝑝,𝐱,Ω)=Ω𝑊𝑝(𝐲)𝑊(𝐲,𝐱)𝑑Ω,𝑛(𝑝,𝐱,Ω)=Ω𝑝(𝐲)𝑊𝑛(𝐲,𝐱)𝑑Ω.(25)

Now with taking into account (23) и (25) integral representations, (22) и (23) can be written in the form𝑤𝑣(𝐱)=𝑊𝑛𝑚,𝐱,𝜕Ω+Θ𝑛𝜃,𝐱,𝜕Ω𝑀𝑛+,𝐱,𝜕Ω𝑉(𝑤,𝐱,𝜕Ω)+𝑊(𝑝,𝐱,Ω)𝐾𝑘=1𝑞𝑐𝑊𝑐(𝐲,𝐱)𝑀𝑐(𝐲,𝐱)𝑤𝑐𝑘,𝐱R2𝜃𝜕Ω𝑛(𝐱)=𝑊𝑛𝑣𝑛,𝐱,𝜕Ω+Θ𝑛𝑚𝑛,𝐱,𝜕Ω𝑀𝑛𝜃𝑛,𝐱,𝜕Ω𝑉𝑛(𝑤,𝐱,𝜕Ω)+𝑊𝑛+(𝑝,𝐱,Ω)𝐾𝑘=1𝑞𝑐𝑊𝑐𝑛(𝐲,𝐱)𝑀𝑐𝑛(𝐲,𝐱)𝜕𝑛𝑤𝑐𝑘,𝐱R2𝜕Ω.(26)

4. Fundamental Solutions

Fundamental solution 𝑊(𝐲,𝐱) for static plate bending is the solution of the differential equation 𝐷ΔΔ(𝐲)𝑊(𝐲,𝐱)=𝛿(𝐲,𝐱),𝐱R2,(27) where ΔΔ(𝐲) is the biharmonic differential operator with respect to 𝐲(𝑦1,𝑦2) defined in (2).

Solution of (27) has the form 1𝑊(𝐲,𝐱)=𝑟8𝜋𝐷2ln𝑟,(28) where 𝑟=(𝑦1𝑥1)2+(𝑦2𝑥2)2 is the distance between points 𝐲 and 𝐱

Fundamental solutions Θ(𝐲,𝐱), 𝑀(𝐲,𝐱) and 𝑉(𝐲,𝐱) from the integral representation (22) can be calculated applying differential operators (3) to 𝑊(𝐲,𝐱) with taking into account (6) and direction of the unit vector 𝐧(𝐲)(see Figure 5).

The fundamental solutions Θ(𝐲,𝐱), 𝑀(𝐲,𝐱), 𝑉(𝐲,𝐱) have the form𝑟Θ(𝐲,𝐱)=8𝜋𝐷𝜕𝑟1𝜕𝑛(1+2ln𝑟),𝑀(𝐲,𝐱)=2𝜕8𝜋2𝑟𝜕𝑛2𝜕+𝜈2𝑟𝜕𝑡2,𝑉𝑟+(1+𝜈)(1+2ln𝑟)(𝐲,𝐱)=4𝜋𝜕𝑟𝜕𝜕𝑛5𝜈22𝑟𝜕𝑛2+𝜕(2𝜈)2𝑟𝜕𝑡2.(29)

Fundamental solutions 𝑊𝑛(𝐲,𝐱), Θ𝑛(𝐲,𝐱), 𝑀𝑛(𝐲,𝐱), and 𝑉𝑛(𝐲,𝐱), from the integral representation (23) can be calculated applying operator of the normal derivative 𝜕𝑛 with respect to 𝐱 to 𝑊(𝐲,𝐱), Θ(𝐲,𝐱), 𝑀(𝐲,𝐱), and 𝑉(𝐲,𝐱) from (28) and (29) with taking into direction of the unit vector 𝐧(𝐱) (Figure 5).

Finally the fundamental solutions 𝑊𝑛(𝐲,𝐱), Θ𝑛(𝐲,𝐱), 𝑀𝑛(𝐲,𝐱), and 𝑉𝑛(𝐲,𝐱) can be represented in the form𝑊𝑛(𝑦𝐲,𝐱)=2(Θ8𝜋𝐷1+2ln𝑟),𝑛1(𝐲,𝐱)=8𝜋𝐷2𝑦2𝑟𝜕𝑟,𝑀𝜕𝑛(1+2ln𝑟)cos𝛾𝑛1(𝐲,𝐱)=4𝜋𝑟22𝑟𝜕𝑟𝜕𝑛cos𝛾𝜈𝜕𝑟𝜕𝑡sin𝛾𝑦2(1+𝜈)+2𝑦2𝜕2𝑟𝜕𝑛2𝜕+𝜈2𝑟𝜕𝑡2,𝑉𝑛1(𝐲,𝐱)=4𝜋𝑟2(𝜈5)𝑟cos𝛾+2𝑦2𝜕𝑟𝜕𝜕𝑛sin𝛾×6𝑟2𝑟𝜕𝑛2cos𝛾+8𝑦2𝜕𝑟×𝜕𝜕𝑛2𝑟𝜕𝑛2𝜕+(2𝜈)2𝑟𝜕𝑡2+2𝑟(2𝜈)+𝜕𝑟𝜕𝑡𝜕𝑟𝜕𝑡cos𝛾2𝜕𝑟.𝜕𝑛sin𝛾(30) Here 𝛾 is the angle between normal vector 𝐧(𝐱) and axis 𝑦2, counted from 𝐧(𝐱) clockwise.

Analysis of (28)–(30) shows that with 𝐱𝐲𝑊(𝐲,𝐱)𝑟2ln𝑟,Θ(𝐲,𝐱)𝑟ln𝑟,𝑀(𝐲,𝐱)ln𝑟,𝑉(𝐲,𝐱)𝑟1,𝑊𝑛(𝐲,𝐱)𝑟ln𝑟,Θ𝑛𝑀(𝐲,𝐱)ln𝑟,𝑛(𝐲,𝐱)𝑟1,𝑉𝑛(𝐲,𝐱)𝑟2.(31) From (31) it follows that for 𝐱𝐲 kernels 𝑊(𝐲,𝐱), Θ(𝐲,𝐱) and 𝑊𝑛(𝐲,𝐱) are continuous, kernels 𝑀(𝐲,𝐱)иΘ𝑛(𝐲,𝐱) are weakly singular, kernels 𝑉(𝐲,𝐱)и𝑀𝑛(𝐲,𝐱) are strongly singular, and kernel 𝑉𝑛(𝐲,𝐱) is hypersingular. Potentials with continuous kernels 𝑊(𝑣𝑛,𝐱,𝜕Ω), Θ(𝑚𝑛,𝐱,𝜕Ω), and 𝑊𝑛(𝑣𝑛,𝐱,𝜕Ω) can be considered in usual (Riemann or Lebegue) sense, potentials with weakly singular kernels 𝑀(𝜃𝑛,𝐱,𝜕Ω) and Θ𝑛(𝑚𝑛,𝐱,𝜕Ω) can be considered as improper, potentials with strongly singular kernels 𝑉(𝑤,𝐱,𝜕Ω), and 𝑀𝑛(𝜃𝑛,𝐱,𝜕Ω) have to be considered in the sense of Cauchy principal values and potential with hypersingular kernel in sense of Hadamard’s finite part 𝑉𝑛(𝑤,𝐱,𝜕Ω) (see [2, 5, 6, 19, 21]).

For 𝐱𝜕Ω above-mentioned singularities may arise in the boundary potentials (24) and therefore not all of them cross the boundary continuously. Boundary properties of these potentials have been well studied in (see [1, 3], etc). Therefore, only final results will be discussed here. These are expressed by the equations 𝑊𝑣𝑛,𝐱,𝜕Ω±𝑣=𝑊𝑛,𝐱,𝜕Ω0,Θ𝑚𝑛,𝐱,𝜕Ω±𝑚=Θ𝑛,𝐱,𝜕Ω0,𝑀𝜃𝑛,𝐱,𝜕Ω±𝜃=𝑀𝑛,𝐱,𝜕Ω0,𝑉(𝑤,𝐱,𝜕Ω)±1=2𝑤(𝐱)+𝑉(𝑤,𝐱,𝜕Ω)0,𝑊𝑛𝑣𝑛,𝐱,𝜕Ω±=𝑊𝑛𝑣𝑛,𝐱,𝜕Ω0,Θ𝑛𝑚𝑛,𝐱,𝜕Ω±=Θ𝑛𝑚𝑛,𝐱,𝜕Ω0,𝑀𝑛𝜃𝑛,𝐱,𝜕Ω±1=±2𝜃𝑛+𝑀𝑛𝜃𝑛,𝐱,𝜕Ω0𝑉𝑛(𝑤,𝐱,𝜕Ω)±=𝑉(𝑤,𝐱,𝜕Ω)0.(32) Here symbols “” and “” denote that two equalities, one with upper and the other with lower signs, are considered. The sign “0” indicates that the direct value of the corresponding potentials on the boundary 𝜕Ω should be taken.

5. Boundary Integral Equations for Boundary-Value Problems

In order to get integral representations for the deflection of the plate 𝑤(𝐱) and slope of the element perpendicular to the middle surface of the plate 𝜃𝑛(𝐱) on the boundary 𝜕Ω limit transition 𝐱𝜕Ω has to be done. With taking into account boundary properties of the elastic potentials (32) on smooth parts of the boundary the following equations can be obtained:12𝑣𝑤(𝐱)=𝑊𝑛𝑚,𝐱,𝜕Ω+Θ𝑛𝜃,𝐱,𝜕Ω𝑀𝑛1,𝐱,𝜕Ω𝑉(𝑤,𝐱,𝜕Ω)+𝑊(𝑝,𝐱,Ω),𝐱𝜕Ω,2𝜃𝑛(𝐱)=𝑊𝑛𝑣𝑛,𝐱,𝜕Ω+Θ𝑛𝑚𝑛,𝐱,𝜕Ω𝑀𝑛𝜃𝑛,𝐱,𝜕Ω𝑉𝑛(𝑤,𝐱,𝜕Ω)+𝑊𝑛(𝑝,𝐱,Ω),𝐱𝜕Ω.(33) For simplicity in this section we did not take into account corner point.

For any boundary-value problem (2), (9)–(11) can be constructed corresponding BIE based on boundary integral representations (33). Several examples of the BIE that correspond to main boundary-value problems for plates are presented bellow.

Boundary-Value Problem (2), (9)

One has𝑊𝑣𝑛𝑚,𝐱,𝜕Ω+Θ𝑛𝑊,𝐱,𝜕Ω+𝑊(𝑝,𝐱,Ω)=0,𝐱𝜕Ω𝑛𝑣𝑛,𝐱,𝜕Ω+Θ𝑛𝑚𝑛,𝐱,𝜕Ω+𝑊𝑛(𝑝,𝐱,Ω)=0,𝐱𝜕Ω.(34)

This system of the BIE is the system of Fredholm integral equations of the first kind with smooth kernels. Such integral equations correspond the so-called ill-posed problems [19]. Different kind of instability can occur during their numerical solution. For their correct solution special regularization technique has to be applied [19].

Boundary Value Problem (2), (10)

One has𝑊𝑣𝑛𝜃,𝐱,𝜕Ω+𝑀𝑛𝑊,𝐱,𝜕Ω+𝑊(𝑝,𝐱,Ω)=0,𝐱𝜕Ω𝑛𝑣𝑛1,𝐱,𝜕Ω2𝜃𝑛(𝐱)𝑀𝑛𝜃𝑛,𝐱,𝜕Ω+𝑊𝑛(𝑝,𝐱,Ω)=0,𝐱𝜕Ω.(35)

This system of the BIE consists of the Fredholm integral equation of the first kind with smooth kernels and singular integral equation. Specific features of the Fredholm integral equation solution have been considered above. Singular integral equations can be solved using technique developed in [14] of the divergent integrals regularization. In process of the system of integral equations (35) solution specific features of both equations have to be taken into account.

Boundary Value Problem (2), (10)

One has12𝜃𝑤(𝐱)+𝑉(𝑤,𝐱,𝜕Ω)+𝑀𝑛,𝐱,𝜕Ω=𝑊(𝑝,𝐱,Ω),𝐱𝜕Ω𝑐12𝜃𝑛(𝐱)+𝑀𝑛𝜃𝑛,𝐱,𝜕Ω+𝑉𝑛(𝑤,𝐱,𝜕Ω)=𝑊𝑛(𝑝,𝐱,Ω),𝐱𝜕Ω𝑐.(36)

Integral operations in this system of integral equations contain kernels with different singularities. Integrals with logarithmic singularity ln𝑟 are weakly singular, they can be consider as improper. Integrals with singularity 𝑟1 are strongly singular and have to be considered in the sense of Cauchy principle value. Integrals with singularity 𝑟2 are hypersingular and have to be considered in the sense of Hadamard’s finite part. For more information about divergent integrals and their application in BIE one can refer to [5, 6, 8, 1316] and references here.

In each specific case corresponding system of BIE can be easily constructed based on (34)–(36).

6. BEM and Approximation of the Region and Functions

The BEM can be treated as the approximate method for the BIE solution, which includes approximation of the functions and the domain where they are defined by discrete finite dimensional model. Let us construct discrete model the plate boundary 𝜕Ω𝑅. We fix in the 𝜕Ω finite number of points 𝐱𝑔(𝑔=1,,𝐺). These points are refered to as global nodes points 𝑉(𝑔)={𝑥𝑔𝜕Ω𝑔=1,,𝐺}.

We shall divide the boundary 𝜕Ω into finite number of subdomains 𝜕Ω𝑛(𝑛=1,,𝑁), such that they satisfy the following conditions: 𝜕Ω𝑚𝜕Ω𝑘=,(ecли𝑚𝑘),𝜕Ω=𝑁𝑘=1𝜕Ω𝑘.(37)

On each FE we introduce a local coordinate system 𝜉. We designate the nodal points 𝐱𝑞𝑉𝑛 in the local system of coordinates by 𝜉𝑞. They are coordinates of the nodal points in the local coordinate system. Local and global coordinate are related in the following way:𝐱𝑞=𝑁𝑛=1Λ𝑛𝜉𝑞𝑛.(38) Functions Λ𝑛 depend on position of the nodal points in the BE. They join individual BE together in the BE model. Borders of the BEs and position of the nodal points should be such that, after joining together, separate elements form discrete model of the plate boundary 𝜕Ω.

Having constructed BE model of the area Ω, we can consider approximation of the function 𝑓(𝐱) that belong to some functional space. The BE model of the boundary 𝜕Ω is the domain of function which should be approximated. We denote function 𝑓(𝐱) on the BE Ω𝑛 by 𝑓𝑛(𝐱). Then𝑓(𝐱)=𝑁𝑛=1𝑓𝑛(𝐱).(39) On each BE the local functions 𝑓𝑛(𝐱) may be represented in the form 𝑓𝑛(𝐱)𝑄𝑞=1𝑓𝑛(𝐱𝑞)𝜑𝑛𝑞(𝜉),(40) where 𝜑𝑛𝑞(𝜉) are interpolation polynomials or shape functions on the BE with number 𝑛. In nodal point with coordinates 𝐱𝑞 they are equal to 1 and in other nodal points are equal to zero. Taking into account (39) and (40) global approximation of the function 𝑓(𝐱) looks like𝑓(𝐱)𝑁𝑄𝑛=1𝑞=1𝑓𝑛𝐱𝑞𝜑𝑛𝑞(𝜉).(41) If the nodal point 𝑞 belongs to several FEs it is considered in these sums only once.

In general case BEs can be of different shape and size with different shape functions defined on them. The simplest are linear BE with piecewise constant shape functions as it is shown in Figure 6.

For linear BE and piecewise constant shape functions corresponding integrals can be calculated analytically. It can significantly simplify calculations, reduce time, and increase accuracy and stability of the calculation process. Advantages of application curvilinear BE and high-order shape functions consist in more accurate approximation of the boundary and functions, but it leads to complication of calculations. Some time that circumstance can devalue the above-mentioned advantages. For more information regarding advantages and disadvantages of different BEs and function approximation refer to [1, 2]. Taking into account all the above we use here the linear BE and the piecewise constant shape functions.

Let us divide boundary 𝜕Ω into 𝑁 curvilinear BE 𝜕Ω𝑛, which satisfy conditions (37). Replacing curvilinear BE by linear ones as it is shown in Figure 6, we obtain discrete approximation of the plate boundary by linear segments of the length 2Δ𝑛. In accordance with model that has been used here 𝑤(𝐱),𝜃𝑛(𝐱), 𝑚𝑛(𝐱)и𝑣𝑛(𝐱) are constant on each BE.

In order to calculate integrals over domain in (25) we divide domain Ω into 𝐿 triangular elements as it is shown in Figure 7.

Therefore, finite dimensional system of the BEM equations has the form 12𝑤𝐱𝑚=𝑁𝑘=1𝑣𝑛𝐲𝑘𝑊𝐱𝑚𝐲𝑘+𝑚𝑛𝐲𝑘Θ𝐱𝑚𝐲𝑘𝜃𝑛𝐲𝑘𝑀𝐱𝑚𝐲𝑘𝐲𝑤𝑘𝑉𝐱𝑚𝐲𝑘+𝑑𝑆𝐿𝑙=1𝑊𝐱𝑚,𝐲𝑙+,𝑝𝐾𝑞𝑐𝑤𝑐,𝐱𝜕Ω𝑛,1(42)2𝜃𝑛𝐱𝑚=𝑁𝑘=1𝑣𝑛𝐲𝑘𝑊𝑛𝐱𝑚𝐲𝑘+𝑚𝑛𝐲𝑘Θ𝑛𝐱𝑚𝐲𝑘𝜃𝑛𝐲𝑘𝑀𝑛𝐱𝑚𝐲𝑘𝐲𝑤𝑘𝑉𝑛𝐱𝑚𝐲𝑘+𝑑𝑆𝐿𝑙=1𝑊𝑛𝐱𝑚,𝐲𝑙+,𝑝𝐾𝑞𝑐𝜕𝑛𝑤𝑐,𝐱𝜕Ω𝑛,(43) where𝑊𝐱𝑚,𝐲𝑘=𝜕Ω𝑘𝑊𝐲,𝐱𝑚Θ𝐱𝑑𝑆,𝑚,𝐲𝑘=𝜕Ω𝑘Θ𝐲,𝐱𝑚Θ𝐱𝑑𝑆,𝑚,𝐲𝑘=𝜕Ω𝑘Θ𝐲,𝐱𝑚Θ𝐱𝑑𝑆,𝑚,𝐲𝑘=𝜕Ω𝑘Θ𝐲,𝐱𝑚𝑊𝑑𝑆,𝑛𝐱𝑚,𝐲𝑘=𝜕Ω𝑘𝑊𝑛𝐲,𝐱𝑚Θ𝑑𝑆,𝑛𝐱𝑚,𝐲𝑘=𝜕Ω𝑘Θ𝑛𝐲,𝐱𝑚𝑑𝑆,(44)𝑀𝑛𝐱𝑚,𝐲𝑘=𝜕Ω𝑘𝑀𝑛𝐲,𝐱𝑚𝑉𝑑𝑆,𝑛𝐱𝑚,𝐲𝑘=𝜕Ω𝑘𝑉𝑛𝐲,𝐱𝑚𝑊𝐱𝑑𝑆,𝑚,𝐲𝑙=,𝑝Ω𝑙𝑝(𝐲)𝑊𝐲,𝐱𝑚𝑊𝑑Ω,𝑛𝐱𝑚,𝐲𝑙=,𝑝Ω𝑙𝑝(𝐲)𝑊𝑛𝐲,𝐱𝑚𝑑Ω.(45)

System of (42) consist of 2𝑁+𝐾 equations, where 2𝑁 corresponds to number of BEs and 𝐾 corresponds to number corner points. For convenience it can be presented in the matrix form 12||||𝑤𝜃𝑛||||=|||||[𝑊Θ]𝑊][𝑛Θ𝑛||||||||||𝑣𝑛𝑚𝑛||||||||||[𝑀𝑉]𝑀][𝑛𝑉𝑛|||||||||𝜃𝑛[𝑤]||||+|||||[𝑊]𝑊𝑛|||||||||[𝑝][𝑝]||||(46) or [𝐴][𝑃].{𝑋}=(47) Here [𝐴] is the matrix of the dimension (2𝑁+𝐾)(2𝑁+𝐾) that depends on boundary conditions, elements of which are functions defined by (44), 𝑋 is column vector of unknown functions of the dimension 2𝑁+𝐾, 𝑃 and is column vector of known functions of the dimension 2𝑁+𝐾 which depend on boundary conditions, external load, and discretization of the domain Ω and its boundary 𝜕Ω.

For any specific external load and boundary conditions from system (46) can be constructed system (47). Matrix 𝐴 and vector 𝑃 can be calculated using (44).

7. Calculation Integrals over Boundary and Domain Elements

There are two approaches to calculation of integrals (44) over the BE. The first one explores numerical calculations using quadrature formulas. The second one consists in analytical integrations of corresponding integrals.

Numerically integrals over boundary elements are usually calculated using the Gaussian quadrature formulas [20, 21]𝑏𝑎𝑝(𝑦)𝑓(𝑦)𝑑𝑥𝑄𝑞=1𝜔𝑞𝑓𝑦𝑞+𝑅(𝑓),(48) where 𝑄 is the number of nodes, [𝑎,𝑏] is any finite or infinite segment of the real line, 𝑦𝑞 are coordinates of the nodes, 𝜔𝑞 is the weighting factor, 𝑝(𝑦) is the weighting function, and 𝑅(𝑓) is the remainder of the quadrature. It is necessary that product of functions 𝑓(𝑦) and 𝑝(𝑦) is integrable on [𝑎,𝑏]. The coordinates 𝑦𝑞 of nodes and the weighting factor 𝜔𝑞 can be found in [20, 21].

The Gaussian quadrature formulas can be effectively applied for calculation of the integrals without singularities, in the case of integrals (44) only for 𝑥𝑚𝜕Ω𝑘. In the case of singular boundary elements, when 𝑥𝑚𝜕Ω𝑘, formulas (48) cannot be applied. In this case special regularization technique has to be applied.

Let us introduce the system of coordinates that related to the BE with number 𝑚; axis 𝑥1 coincides with that BE and axis 𝑥2 is perpendicular to it, as it is shown in Figure 6. Calculation of the regular and divergent integral will be done in that system of coordinates.

7.1. Calculation of the Divergent Integrals

In the above mentioned system of coordinates, the local coordinates of the points 𝐱𝑘 and 𝐲𝑚 have the form𝑥𝑚1=0,𝑥𝑚2=0,𝑦𝑚1=𝑡,𝑦𝑚2=0,(49) where 𝑡[Δ𝑚,Δ𝑚] is a local coordinate.

Therefore, fundamental solutions on the singular boundary element have the form [17]𝑊𝐱𝑚=1,𝐲𝑡8𝜋𝐷2𝑀𝐱ln|𝑡|,𝑚1,𝐲=[],Θ8𝜋2𝜈+(1+𝜈)(1+2ln|𝑡|)𝑛𝐱𝑚1,𝐲=𝑉8𝜋𝐷(1+2ln|𝑡|),𝑛𝐱𝑚=,𝐲1+𝜈4𝜋𝑡2,Θ𝐱𝑚𝐱,𝐲=𝑉𝑚,𝐲=𝑊𝑛𝐱𝑚,𝐲=𝑀𝑛𝐱𝑚,𝐲=0.(50)

Integrals in (50) are not complicate and can be calculated analytically using regularization technique developed in [16]. As result we obtain 𝑊𝐱𝑚,𝐲𝑚=Δ3𝑚36𝜋𝐷3lnΔ𝑚,𝑀𝐱1𝑚,𝐲𝑚Δ=𝑚4𝜋𝜈1+2(1+𝜈)lnΔ𝑚,Θ𝑛𝐱𝑚,𝐲𝑚Δ=𝑚4𝐷2lnΔ𝑚,𝑉1𝑛𝐱𝑚,𝐲𝑚=1𝜈2𝜋Δ𝑚.(51)

In (51) it has been taken into account that integral containing kernel 𝑊 is regular and no special consideration is needed. The integrals containing kernels 𝑀 and Θ𝑛 are strongly singular; they have to be considered in the sense of Cauchy principle value. The integral containing kernel 𝑉 is hypersingular; it has to be considered in the sense of Hadamard’ finite part. For more information refer to [13, 14, 22].

7.2. Calculation of the Regular Integrals

In the system of coordinates presented in Figure 6 coordinates of point 𝐱𝑘 and 𝐲 are𝑥𝑚1=0,𝑦1=𝑦𝑘1𝑥+𝑡cos𝛾,𝑚2=0,𝑦2=𝑦𝑘2𝑡sin𝛾,(52) where 𝑦𝑘1, 𝑦𝑘2 are coordinates of the middle of the BE with number 𝑘, 𝛾 is the angle between vectors 𝐧𝑚 and 𝐧𝑘 that are perpendicular to the BEs with numbers 𝑚 and 𝑘, respectively counted in clockwise direction. Fundamental solutions in this case are regular and are presented by (28)–(30). After analytical integration they can be presented in the following form [17]:𝑊𝐱𝑚,𝐲𝑘=1𝑟8𝜋𝐷2𝐼1+2𝐵𝐼2+𝐼3,Θ𝐱𝑚,𝐲𝑘=𝐶𝐼4𝜋𝐷1+Δ𝑘,𝑀𝐱𝑚,𝐲𝑘1=𝐼4𝜋(1+𝜈)1+Δ𝑘+𝐼1,0𝐶2+𝜈𝐵2+2𝐼1,1𝜈𝐵+𝐼1,2𝜈,𝑉𝐱𝑚,𝐲𝑘𝐶=𝐼4𝜋1,0(5𝜈)2𝐼2,0𝐶2+𝜈𝐵2(2𝜈)2(2𝜈)2𝐼2,1𝜈𝐵+𝐼2,2,𝑊𝐱𝑚,𝐲𝑘1=𝑦4𝜋𝐷2𝐼1+Δ𝑘𝐼2,Θsin𝛾𝑛𝐱𝑚,𝐲𝑘1=×𝐼4𝜋𝐷1+Δ𝑘𝐼cos𝛾+𝐶1,0𝑦2𝐼1,1;sin𝛾(53)𝑀𝑛𝐱𝑚,𝐲𝑘1=𝐼4𝜋1,02𝐵sin𝛾2𝐶cos𝛾(1+𝜈)𝑦2+𝐼1,1(1+3𝜈cos𝛾)+2𝐼2,0𝑦2𝐶2+𝐵2𝜈+2𝐼2,0𝑦2𝐶2+𝐵2𝜈+2𝐼2,1×2𝐵𝜈𝑦2𝐶2+𝐵2𝜈sin𝛾+2𝐼2,2𝜈𝑦22𝐵sin𝛾,2𝐼2,3,𝑉𝜈sin𝛾𝑛𝐱𝑚,𝐲𝑘1=𝐼4𝜋1,0(𝜈5)cos𝛾+2𝐼2,0𝐶(5𝜈)𝑦2+𝐵(2𝜈)×(𝐵cos𝛾2𝐶sin𝛾)+3𝐶2cos𝛾+2𝐼2,1(𝐶(5𝜈)sin𝛾+2(2𝜈)×(𝐵cos𝛾𝐶sin𝛾))+2𝐼2,2(2𝜈)×cos𝛾8𝐼3,0𝐶𝑦2𝐶2+(2𝜈)𝐵28𝐼3,1𝐶𝐶2×sin𝛾+𝐵(2𝜈)2𝑦2𝐵sin𝛾8𝐼3,2×𝑦(2𝜈)22𝐵sin𝛾×𝐶+8𝐼3,3.𝐶(2𝜈)sin𝛾(54) Here 𝐵=𝑦1cos𝛾𝑦2sin𝛾,𝐷=arctg(Δ𝑘+𝐵)/𝐶arctg(Δ𝑘+𝐵)/𝐶,𝐶=𝑦1sin𝛾+𝑦2cos𝛾,𝐹𝐹=1𝐹2||Δ=ln2𝑘+2𝐵Δ𝑘+𝑟2||||Δ2𝑘2𝐵Δ𝑘+𝑟2||,||||=𝑦𝑟=𝐲𝐱21+𝑦221/2,𝐼𝑛,𝑚=Δ𝑘Δ𝑘𝑡𝑚𝑑𝑡𝑡2+2𝐵𝑡+𝑟2𝑛𝐼,𝑛=1,2,3,𝑚=1,2,3,𝑛=Δ𝑘Δ𝑘𝑡𝑛ln𝑟𝑑𝑡,𝑛=0,1,2.(55)

7.3. Calculation of the Domain Integrals

In the BIE can appear also domain integrals (25), which are used to take into account action of the load distributed over domain on the plate. They are regular and can be calculated numerically using the Gaussian quadratures [20, 21].

In order to calculate integrals over domain Ω we divide it into 𝐿 triangular elements as it is shown in Figure 7. Integrals in (45) over triangle can be calculated in the form 𝑊𝐱𝑚,𝐲𝑙=,𝑝𝑄𝑞=1𝑝𝐲𝑞𝑊𝐱𝑚,𝐲𝑞𝜔𝑞,𝑊𝑛𝐱𝑚,𝐲𝑙=,𝑝𝑄𝑞=1𝑝𝐲𝑞𝑊𝑛𝐱𝑚,𝐲𝑞𝜔𝑞,𝐱𝑚Ω,(56) where 𝑦𝑞 are coordinates of nodes on the triangular element, 𝜔𝑞 is the weighting factor. Coordinates 𝑦𝑞 of nodes and weighting factor 𝜔𝑞 for triangle can be found in [20, 21].

8. Numerical Examples

We have considered here some benchmark examples that correspond to bending of the thin plates of different shape. In all examples the ratio 2/𝑅=0.1 has been used.

8.1. Calculation of the Circular Plate

First let us consider circular simply supported over-the-contour plate that subjected to action the uniformly distributed over the domain Ω load. For this case analytical solution has the form [23] 𝑤(𝜌)=𝑃𝜌2𝜌16𝐷(1+𝜈)2𝑅2,(57) where 𝜌 is the radial polar coordinate, 𝑅 is the radius of the plate.

Expressions for 𝜃𝑛, 𝑚𝑛, 𝑣𝑛 can be obtained from (57) applying differential operators (3). Dependence of the accuracy of the BEM on number 𝑁 of the BE has been studied here. In all occurrences for approximation of the domain Ω it was divided by 288 triangular elements. Results of calculations of the unknown boundary data are presented on Figure 8, where curve 1 corresponds to analytical solution, curve 2 corresponds to Δ𝜃𝑛=𝜃appr𝑛/𝜃extact𝑛, and curve 3 corresponds to Δ𝑣𝑛=𝑣appr𝑛/𝑣extact𝑛, respectively.

Analysis of the data presented in Figure 8 shows that generalized shear force 𝑣𝑛 is calculated more accurately. For example, calculation accuracy to within 5% for 𝑣𝑛 is reached for 𝑁=13, whereas the same accuracy for the slope of the element perpendicular to the middle surface 𝜃n is reached for 𝑁=24.

In order to compare traditional and proposed here analytical approach, evaluation of integrals (44) using the Gaussian quadratures (48) and presented here equations (52) has been done for circular simply supported over-the-contour plate. In both cases divergent integrals have been calculated using (50) and 𝑁=32 rectilinear BE. Time for calculation of the above-mentioned integrals in traditional approach has to be four times longer in order to reach the same accuracy as in proposed here analytical approach.

8.2. Calculation of the Rectangular Plate

We consider here some benchmark examples for rectangular plate with different boundary conditions. In the first example plate was loaded by concentrated force applied to the point in its center. In all other examples the plate was loaded by uniformly distributed over the domain Ω load. Boundary conditions and load are presented in the first column of Tables 16. For all presented here examples solution obtained by BEM was compared with existing analytical solution, which can be presented in the form of Fourier series [23].

Loading patternNumber of BEPointsBoundary conditions Δ 𝑡 Δ 𝑤 m a x
Δ 𝑤 Δ 𝜃 𝑛 Δ 𝑀 𝑛 Δ 𝑉 𝑛


Loading patternNumber of BEPointsBoundary conditions Δ 𝑡 Δ 𝑤 m a x
Δ 𝑤 Δ 𝜃 𝑛 Δ 𝑀 𝑛 Δ 𝑉 𝑛


Loading patternNumber of BEPointsBoundary conditions Δ 𝑡 Δ 𝑤 m a x
Δ 𝑤 Δ 𝜃 𝑛 Δ 𝑀 𝑛 Δ 𝑉 𝑛


Loading patternNumber of BEPointsBoundary conditions Δ 𝑡 Δ 𝑤 m a x
Δ 𝑤 Δ 𝜃 𝑛 Δ 𝑀 𝑛 Δ 𝑉 𝑛


Loading patternNumber of BEPointsBoundary conditions Δ 𝑡 Δ 𝑤 m a x
Δ 𝑤 Δ 𝜃 𝑛 Δ 𝑀 𝑛 Δ 𝑉 𝑛


Loading patternNumber of BEPointsBoundary conditions Δ 𝑡 Δ 𝑤 m a x
Δ 𝑤 Δ 𝜃 𝑛 Δ 𝑀 𝑛 Δ 𝑉 𝑛


Columns 4–7 of Tables 16 show Δ𝑤,Δ𝜃𝑛,Δ𝑚𝑛,Δ𝑣𝑛—the ratio of the unknown boundary conditions for different values of 𝑁 to the values obtained analytically at points 1 and 2 for the schemes presented in the column 1. This data can be used to optimize the number of the BE depending on the required accuracy for different boundary conditions. According to our estimates, the value 𝑁=40 can be recommended in calculations performed for different combinations of boundary conditions and loading schemes. Here, the error in the calculation of the unknown boundary conditions does not exceed 1%. Column 9 of Tables 16 shows values of Δ𝑤max—ratio of the maximum deflection calculated by the BEM to the corresponding values obtained analytically [23]. With number of the BE equal to 𝑁=40, the error is no greater than 1.7% for any variant of boundary conditions and loading presented in Tables 16 except the first one. It is equal to 4.5% for the first case, which can be attributed to the inaccuracy of approximation of the concentrated force.

Integrals in (44) have been calculated numerically using the Gaussian quadratures (48) and analytically using presented here equations (52). An important fact testifying to the expediency of use of the presented here equations for analytical integration (52) for calculation of the coefficients of system (47) is the reduction of the computational time. Column 8 of the Tables 16 shows values of Δ𝑡—the ratio of the time required to construct the system of the linear algebraic equations (47) analytically using (52) and numerically using the Gaussian quadratures. Studies have shown that use of the approach proposed here makes it possible to cut the time required to obtain solution the problem by the factor from 2.5 to 5. This is particularly important when the calculations have to be performed repeatedly, such as in the solution of the dynamic and nonlinear problems. The amount of time saved depends mainly on the number of the BE approximating the boundary, as well as on the type of boundary conditions.

In order to visualize results from Tables 16 also Figures 914 are presented. They complement data from the corresponding tables and present them in more convenient and visual form.

8.3. Calculation of the T-Shaped Plate

The diagram of the T-shaped plate together with boundary conditions is presented in Figure 15. The plate is subjected to action of the uniformly distributed over the domain Ω load of the magnitude 𝑞. Material properties of the plate have been calculated from the equation (𝑞𝑎4)/(𝐸4)=100 and Poisson ratio 𝜐=0.3.

In Figure 16 are presented distributions of the values 𝜎𝑥𝑎2/4𝐸2 and 𝑤/ for the cross-section with coordinate 𝑥2=0 and for number of the Bes𝑁=96. The solid line corresponds to the BEM and dash line to FEM solutions, respectively.

From this data follows that results obtained by the BEM and by the FEM are in a good agreement but time of calculation by the BEM is significantly less.

9. Conclusion

Direct BIEM based on Betty’s theorem is applied here for solution of thin elastic plate bending problems for different boundary conditions and load. Analytical integration of the regular and divergent integrals over the BE is applied, and effective formulas for calculation of coefficients of the system of linear algebraic equations for the BEM have been developed. The main advantage of the proposed approach consists in significant reduction of the calculation time comparison with traditional approach based on Gaussian’s quadratures. Numerical examples demonstrate effectiveness of the proposed here approach. In all presented examples it was demonstrated high accuracy, in good agreement with existing analytical solutions and significant reduction of the time of calculations in comparison with traditional approaches.


The author is very grateful to his former Ph.D. student Dr. Alexander Lukin from Kharkov State University for help in this paper preparation.


  1. M. H. Aliabadi, “Applications in solids and structures,” in The Boundary Element Method, vol. 2, John Wiley &Sons, 2002. View at: Google Scholar
  2. F. Hartmann, Introduction to Boundary Elements. Theory and Applications, Springer, 1989.
  3. D. E. Beskos, Ed., Boundary Element Analysis of Plates and Shells, Springer, Berlin, Germany, 1991.
  4. N. J. Altiero and D. L. Sikarskie, “A boundary integral method applied to plates of arbitrary plan form,” Computers and Structures, vol. 9, pp. 163–168, 1978. View at: Google Scholar
  5. A. Frangi and M. Guiggiani, “Boundary element analysis of Kirchhoff plates with direct evaluation of hypersingular integrals,” International Journal for Numerical Methods in Engineering, vol. 46, pp. 1845–1863, 1999. View at: Google Scholar
  6. A. Frangi and M. Guiggiani, “A direct approach for boundary integral equations with high-order singularities,” International Journal for Numerical Methods in Engineering, vol. 49, pp. 871–898, 2000. View at: Google Scholar
  7. F. Hartmann and R. Zormantel, “The direct boundary element method in plate bending,” International Journal for Numerical Methods in Engineering, vol. 23, pp. 2049–2069, 1986. View at: Google Scholar
  8. B. Knopke, “The hypersingular integral equation for the bending moments mxx, mxy, and myy of the Kirchhoff plate,” Computational Mechanics, vol. 15, pp. 19–30, 1994. View at: Google Scholar
  9. M. Stern, “A general boundary integral formulation for plate bending problems,” International Journal for Solids and Structures, vol. 15, pp. 769–782, 1979. View at: Google Scholar
  10. P. H. Wen, M. H. Aliabadi, and A. Young, “Application of dual reciprocity method to plates and shells, engineering analysis with boundary element,,” International Journal for Numerical Methods in Engineering, vol. 24, pp. 583–590, 2000. View at: Google Scholar
  11. V. V. Zozulya and A. N. Lukin, “Calculation of plates by the method of boundary integral equations under the dynamic loading,” Doclady Akademii Nauk of Ukraine, vol. 8, pp. 64–70, 1997. View at: Google Scholar
  12. V. V. Zozulya and A. N. Lukin, “Calculation of plates by the boundary integral equations method,” Doclady Akademii Nauk of Ukraine, vol. 4, pp. 39–45, 1996. View at: Google Scholar
  13. J. Hadamard, Lectures on Cauchy's Problem in Linear Partial Differential Equations, Yale University Press, New Haven, Conn, USA, 1923.
  14. I. K. Lifanov et al., Hypersingular Integral Equations and Their Applications, Chapman&Hall, 2004.
  15. V. V. Zozulya, “Regularization of the divergent integrals. I. General consideration,” Electronic Journal of Boundary Elements, vol. 4, no. 2, pp. 49–57, 2006. View at: Google Scholar
  16. V. V. Zozulya, “The Regularization of the Divergent Integrals in 2-D Elastostatics,” Electronic Journal of Boundary Elements, vol. 7, no. 2, pp. 50–88, 2009. View at: Google Scholar
  17. V. V. Zozulya and A. N. Lukin, “Calculation of plates by the boundary element method,” Doclady Akademii Nauk of Ukraine, vol. 4, pp. 60–67, 1996. View at: Google Scholar
  18. M. A. Krasnosel'skii et al., Integral Operators in Spaces of Summable Functions, Noordhoff, 1976.
  19. A. N. Tikhonov and V. Y. Arsenin, Solutions of Ill-Posed Problems, Winston, 1977.
  20. P. Kythe and M. Schaferkotter, Handbook of Computational Methods for Integration, CRC, 2005.
  21. V. Krylov, Approximate Calculation of Integrals, Dover, 2006.
  22. V. V. Zozulya and A. N. Lukin, “Design of plates by the boundary element method,” International Applied Mechanics, vol. 33, pp. 241–244, 1997. View at: Google Scholar
  23. S. Timoshenko and S. Woinowski-Krieger, Theory of Plates and Shells, McGraw-Hill, New York, 1959.

Copyright © 2011 V. V. Zozulya. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Related articles

No related content is available yet for this article.
 PDF Download Citation Citation
 Download other formatsMore
 Order printed copiesOrder

Related articles

No related content is available yet for this article.

Article of the Year Award: Outstanding research contributions of 2021, as selected by our Chief Editors. Read the winning articles.