- Email: [email protected]
BEM formulation for von Kármán plates
ARTICLE IN PRESS Engineering Analysis with Boundary Elements 33 (2009) 1223–1230
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
BEM formulation for von Ka´rma´n plates Leandro Waidemam a, Wilson S. Venturini b, a b
´ Federal University of Technology, BR 369, km 0.5, 87301-006-Campo Moura ˜o, Brazil Parana ˜o Carlos School of Engineering, University of Sa ˜o Paulo, Av. Trabalhador Sa ˜o-Carlense 400, 13566-970 Sa ˜o Carlos, Brazil Sa
a r t i c l e in fo
abstract
Article history: Received 1 April 2007 Accepted 24 April 2009 Available online 7 June 2009
This work deals with nonlinear geometric plates in the context of von Ka´rma´n’s theory. The formulation is written such that only the boundary in-plane displacement and deflection integral equations for boundary collocations are required. At internal points, only out-of-plane rotation, curvature and in-plane internal force representations are used. Thus, only integral representations of these values are derived. The nonlinear system of equations is derived by approximating all densities in the domain integrals as single values, which therefore reduces the computational effort needed to evaluate the domain value influences. Hyper-singular equations are avoided by approximating the domain values using only internal nodes. The solution is obtained using a Newton scheme for which a consistent tangent operator was derived. & 2009 Elsevier Ltd. All rights reserved.
Keywords: Bending plates Geometrical nonlinearities
1. Introduction The boundary element method (BEM) has been successfully applied to solve plate bending problems many times. An important characteristic of the boundary methods applied to plate bending is that all boundary values are approximated by the same shape function, which avoids the necessity of using higher order derivatives of the displacement approximation to compute internal forces. Thus, bending and twisting moments and shear forces are precisely evaluated. The method has already proved to be accurate and reliable enough for this kind of application. The first works using the indirect version of the boundary methods applied to Kirchhoff’s plates were provided by Jaswon et al. [1] and Jaswon and Maiti [2]. The direct formulation based on the same simplified plate bending hypothesis was proposed by Bezine [3], Tottenhan [4] and Stern [5]. Van der Weee¨n [6] has derived the boundary element formulation in the context of the Reissner’s theory of plates. Barcellos and Silva proposed the extension to consider the Mindlin hypothesis [7]. More recent contributions in boundary elements to treat shear deformable plates can be found in Aliabadi [8]. The plate bending numerical formulation is a very important subject in engineering because of its many applications to a large number of complex structures, such as aircraft, ships, cars, pressure vessels, and off shore structures, among others. Usually, these complex problems require accurate plate bending models to
Corresponding author.
E-mail addresses: [email protected] (L. Waidemam), [email protected] (W.S. Venturini). 0955-7997/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2009.04.012
take into account physical and geometric nonlinear effects. Regarding the use of plastic models in conjunction with plate bending BEM models, we have to point out the works carried out by Moshaiov and Vorus [9]. Chueiri and Venturini [10] who have proposed a general model that takes into account not only the plastic criterion written in terms of moments, by also the integration of the corrector stress field over the thickness, which is very convenient for the modelling of reinforced concrete behaviour. Providakis [11] was the first to propose a viscoplastic model for plate bending analysis. Several BEM formulations have already been proposed to consider nonlinear geometric effects in the global equilibrium of plate bending problems. One of the first works treating this subject was proposed by Morjaria [12]. Kamiya and Sawaki [13] have proposed a BEM formulation for elastic plates governed by the Berger equation. The first BEM formulation to analyse plate bending problems within the context of von Ka´rma´n hypothesis was proposed by Ye and Liu [14], who used a fictitious load distributed over the domain to model nonlinear effects. The von Ka´rma´n hypothesis was also adopted by Tanaka et al. [15] to develop a more elaborated BEM incremental formulation to deal with finite deflections of thin elastic plates. Wang et al. [16] have worked on von Ka´rma´n plates and introduced the dual reciprocity approach based on global radial functions to approximate the correcting integral term. All works reported above were proposed in the context of thin plates. Several other important works have appeared more recently that demonstrate the efficiency with which BEM formulations deal with shear deformable plates based on the Reissner–Mindlin hypothesis and taking into account geometric and material nonlinearities. Among them, we emphasise the works carried out by Aliabadi and his
ARTICLE IN PRESS 1224
L. Waidemam, W.S. Venturini / Engineering Analysis with Boundary Elements 33 (2009) 1223–1230
co-workers: Wen, et al. [17], Purbolaksono and Aliabadi [18], Dirgantara and Aliabadi [19], and Supriyono and Aliabadi [20]. In this work, we return to the BEM formulation based on von Ka´rma´n’s theory. Emphasis is placed on the accurate evaluation of the domain integrals that are approximated with cells and to the solution technique for which a tangent consistent operator is proposed. Examples of plates with finite deflection are analysed, and the results are compared with other numerical solutions.
þ
Nc X
Rc ðC k Þwc ðS; C k Þ þ
k¼1
Z þ
O
Z
gðpÞw ðS; pÞdO Og
N ij ðpÞw;ij ðpÞw ðS; pÞ dO Z
(3)
P ij ðS; PÞuj ðPÞ dG þ
Z
uij ðS; PÞpj ðPÞ dG Z Z 1 þ uij ðS; pÞbj ðpÞ dO N ðS; pÞw;j ðpÞw;k ðpÞ dO 2 O ijk O
C ij ðSÞuj ðSÞ ¼
G
G
(4)
2. Theoretical background Without a loss of generality, let us consider a single thin plate region O with boundary G over which a distributed load q is applied orthogonally to the middle surface, i.e. in the direction x3, as shown in Fig. 1. This plate region can also be subjected to in-plane forces (directions x1 and x2) that are either distributed over the domain or applied along the boundary. To write the field equations for this plate problem, the following hypothesis can be assumed according to the von Ka´rma´n theory: the strains are assumed to be small and the final deflections are on the order of the plate thickness h. For this problem, the differential equations of the bending and in-plane problems are given by [21] mij ;ij þ N ij w;ij bi w;i þ mi ;i þ q ¼ 0
(1a)
Nij ;j þ bi ¼ 0
(1b)
where w represents the deflections, mij are bending and twisting moments, Nij gives the in-plane internal forces, bi represents the in-plane domain loads applied and mi is the applied moment over the plate domain; the subscripts are in the range i,j ¼ {1, 2}. The internal forces in Eqs. (1a) and (1b) can be written in terms of displacements as follows: mij ¼ Dðndij w;kk þ ð1 nÞw;ij Þ
where ( )* represents the well-known fundamental solutions of both problems; Vn, Mn, w and w,n are the effective shear forces, moments, deflections and rotations orthogonal to the plate boundary applied along the boundary, respectively; Rc and wc represent the corner reactions and deflections, respectively; pj and uj represent, respectively, the in-plane tractions and displacements along the plate boundary; w,j, w,ij and Nij are the rotations, curvatures and in-plane forces, respectively; S, P, p and Ck represent the collocation point and the field points defined along the boundary, inside the domains and at the corners, respectively; C(S) and Cij(S) are the well-known free terms of the plate bending and stretching problems [22]. To complete the required integral equations, one has to obtain the integral representations of rotations, curvatures and internal forces for a given domain collocation s as follows: Z w;i ðsÞ ¼ fV n ;i ðs; PÞwðPÞ M n ;i ðs; PÞw;n ðPÞg dG G
Nc X k¼1
þ
Nc X
Eh ð1 nÞ
2
nðuk ;k þ w;k w;k =2Þdij þ
Z
(2a) 1n ðui ;j þ uj ;i þ w;j w;i Þ 2
w;ij ðsÞ ¼ (2b)
CðSÞwðSÞ þ
Z G
Z ¼ G
ðV n ðS; PÞwðPÞ M n ðS; PÞw;n ðPÞÞ dG þ
Nc X
Rc ðS; C k Þwc ðC k Þ
k¼1
ðV n ðPÞw ðS; PÞ Mn ðPÞw;n ðS; PÞÞ dG
Γ x1 x2
Ω
x3
Fig. 1. General plate domain.
Z G
Nc X
Rc ;ij ðs; C k Þwc ðC k Þ þ
Nc X
Rc ðC k Þwc ;ij ðs; C k Þ þ
k¼1
Z þ
O
Nij ðsÞ ¼
G
fV n ðPÞw;i ðs; PÞ M n ðPÞw;ni ðs; PÞg dG
Z Og
gðpÞw;i ðs; pÞ dO
(5)
fV n ;ij ðs; PÞwðPÞ M n ;ij ðs; PÞw;n ðPÞg dG
k¼1
þ
Z
N jk ðpÞw;jk ðpÞw;i ðs; pÞ dO
þ
where ui are the in-plane displacement components. The problem definition is then completed by assuming the following boundary conditions over G: ui ¼ u% i on G1 (generalised displacements, in-plane displacements, deflections and rotations); and pi ¼ p¯i on G2 (generalised tractions, in-plane tractions normal bending moment and effective shear forces), where G1[G2 ¼ G. To obtain the integral equations of the plate bending and stretching problems considering geometrical nonlinearities within the context previously defined, one can apply the Betti’s reciprocity to the linear parts of the stress and strain fields. The integral representations of both problems are given by [22]:
Rc ðC k Þwc ;i ðs; C k Þ þ
k¼1
O
Nij ¼
Rc ;i ðs; C k Þwc ðC k Þ þ
Z G
fV n ðPÞw;ij ðs; PÞ M n ðPÞw;nij ðs; PÞg dG
Z Og
gðpÞw;ij ðs; pÞ dO
N km ðpÞw;km ðpÞw;ij ðs; pÞ dO
Z Dijk ðs; PÞpk ðpÞ dG Sijk ðs; PÞuk ðPÞ dG G G Z Z 1 þ Dijk ðs; pÞbk ðpÞ dO T ðs; pÞw;k ðpÞw;l ðpÞ dO 2 O ijkl O Gh ð2w;i ðsÞw;j ðsÞ þ w;k ðsÞw;k ðsÞdij Þ þ 8ð1 u Þ
(6)
Z
(7)
where v¯ ¼ v/(1+v) is used to simulate the plane stress conditions; the higher order tensors are given in the appendix. The free term appearing in Eq. (7) is due to the differentiation of a singular integral. It has been derived for 3D, plane stress and plane strain problems by Bui [23] and extended for several other integral equations containing kernels with the same singularity.
Ωg 3. Algebraic equations Eqs. (3)–(7) were obtained from von Ka´rma´n’s field equations; therefore, they are written in terms of total values. To reduce the number of domain integrals in this section, it is convenient to replace the densities given by the product of the two values,
ARTICLE IN PRESS L. Waidemam, W.S. Venturini / Engineering Analysis with Boundary Elements 33 (2009) 1223–1230
which appear in the domain integrals, by single equivalent values. Then, the new equivalent values will be approximated by using node values and the shaped functions. Thus, the domain value Nij(p)w,ij(p) at a field point p that appears in the bending problem Eqs. (3), (5) and (6) will be replaced by an equivalent scalar value T(p). Similarly, in Eqs. (4) and (7) of the stretching problem, we replace the domain density w,j(p)w,k(p) by an equivalent tensor Wjk(p). As is typical when solving a nonlinear problem using an incremental process, it is convenient to write the integral representations derived in the previous section into their corresponding rate integral representations. Bearing in mind that we have already replaced the densities in the bending problem equations by the equivalent scalar T(p) and in the in-plane problem equations by the equivalent tensor Wjk(p), we now have to ˙ jk(p), replace these values by the corresponding rates T˙(p) and W respectively. All the other boundary and domain values appearing in all integral equations derived in the previous sections are also replaced by their rates. The starting point of the numerical formulation for the formulated boundary value problem is its incremental representation. The governing equations and the derived integral representations have to be written in increments. To find their appropriate incremental forms, we assume that Dt ¼ tn+1tn is a typical time-step in the time discretization. Then, any rate quantity x˙ integrated along this time interval Dt becomes Dx ¼ xn+1xn, which will replace x˙ in all governing equations derived from the integral representations. Thus, the rate integral representations (3)–(7) are transformed into incremental representations. For instance, Eq. (3) will be first written in rates and then in increments as follows: CðSÞDwðSÞ þ
Z G
þ
Nc X
ðV n ðS; PÞDwðPÞ Mn ðS; PÞDw;n ðPÞÞ dG
Rc ðS; C k ÞDwc ðC k Þ
k¼1
Z ¼
G
þ
ðDV n ðPÞw ðS; PÞ DM n ðPÞw;n ðS; PÞÞ dG
Nc X
DRc ðC k Þwc ðS; C k Þ þ
k¼1
Z
DTðpÞw ðS; pÞ dO
þ
Z
DgðpÞw ðS; pÞ dO
Og
(8)
O
The same modification has to be performed on the other integral representations, Eqs. (4)–(7). After obtaining the finite step integral equations, one has to transform them into algebraic equations by approximating the boundary and domain incremental values. For this purpose, the plate boundary G is discretized into a series of elements, Gs, along which generalised displacements and tractions are approximated by using continuous and discontinuous linear boundary elements. For the present plate stretching–bending problem, discontinuities of boundary values are always expected, particularly at corners and traction jumps. The boundary value discontinuities are always introduced by defining the collocation points placed inside the elements. Typically, for the plate stretching and bending boundary element formulations, we need two equations for each problem for each boundary node used to discretize the boundary. Thus, two algebraic equations for each node are referred to the stretching problems and derived from the incremental form of Eq. (4). Another two algebraic equations for the plate bending problem are derived from the deflection Eq. (8). Since we are avoiding writing the slope equation for this purpose, we have to write two independent deflection integral representations for
1225
each node; therefore, two collocations for each boundary node have to be defined: one on the boundary, either defined at the node or inside the element, and the second out of the domain. Moreover, the corner reactions can also be properly accounted for in the effective shear forces of the adjacent boundary elements, avoiding an additional unknown corner value. This scheme has proved to lead to stable and accurate results. As previously demonstrated, the required outside collocations have to be defined at points very near the boundary to guarantee the accuracy and stability of the solution [24,25]. For this work, they have been defined within the range of 0 to half of the adjacent element average length. It is important to stress that corner double nodes are always used to guarantee the discontinuities of all boundary values. Thus, the boundary values that refer to nodes defined at both elements adjacent to a corner are always independent of each other. Moreover, extra values of the deflection and the concentrated reaction can be associated with the corner. If these values are assumed to be independent, one extra equation for this corner has to be written. In this work, we prefer to replace the corner deflection by the average deflection values of the two double nodes, while the concentrated reaction is properly distributed along the adjacent elements and added to the effective shear force [26]. Accurately performing the integrals over the elements requires special care, particularly when carrying out quasi-singular integrals. To assure the accuracy of these integrals, an appropriate sub-element scheme was employed. This sub-element scheme has already been successfully used to compute quasi-singular integrals that result from the displacement and stress equations for linear and nonlinear stretching problems [27–29]. The accuracy of this scheme has also been confirmed when applied to plate bending problems [30,31]. In these works, the accuracy of the adopted scheme that performs numerically quasi-singular integrals was certified by comparing the results with the values obtained by using analytical expressions and also by verifying the rigid body movements. To obtain the algebraic forms of the Eqs. (4)–(8) the domain is also discretized into continuous and discontinuous linear cells over which the relevant domain values are approximated, wherever required. Therefore, the domain integrals are all approximated over the cells using the appropriate shape functions and nodal values. For these approximations, only domain values are used to approach the domain fields; therefore, discontinuous cells are always required for cells adjacent to the boundary (Fig. 2). Moreover, boundary and domain discretizations may be completely independent of each other. Defining cells with sides larger than the corresponding adjacent boundary elements is recommended to avoid choosing internal nodes too close to the boundary. To describe the scheme used to compute the integrals over a given cell Om, let us take the last domain integral in Eq. (8). Assuming the density DT(p) is approximated over the given cell Om using the shape function Fk(p) and the nodal values DTwk, one can write Z
DTðpÞw ðs; PÞ dO ¼
Om
Z Om
Fk ðpÞw ðs; PÞ dODT kw
(9)
For the linear approximation case, the shape function Fk at a given field point p, written in terms of a global Cartesian system of coordinates, is
Fm ðpÞ ¼
1 ½am þ mm xðpÞ þ nm yðpÞ 2A
(10)
where A is the cell area and the constants am, mm and nm are the well-known coefficients of the linear approximation [22].
ARTICLE IN PRESS 1226
L. Waidemam, W.S. Venturini / Engineering Analysis with Boundary Elements 33 (2009) 1223–1230
Fig. 2. Linear domain approximation: (a) over continuous cells; (b) over discontinuous cells.
represents the increment of the product of the in-plane forces by curvatures. Although the integrals along the cell sides in Eq. (12) are not singular, one can use the same proposed sub-element scheme to compute the quasi-singular boundary elements to improve the computed vales. After performing all boundary and cell integration and storing the coefficients into the corresponding matrices, Eq. (8) becomes
Γ ξ2
Ω
2
w
w w b b w Hw b DU ¼ Gb DP þ S b DT N þ DBb
p
b
where DU and DP contain all boundary displacement and traction nodal values of the bending problem, plus the additional corner displacements and reactions, DTNw represents the condensed summation of the in-plane forces multiplied by the curvatures, Sbw is the corresponding matrix obtained by integrating all domain cells and DBbw gives the domain load effects. Similarly, one can also transform the in-plane displacement integral representation, Eq. (4), which has already modified by replacing the boundary and internal values by their increments, into its algebraic form. We first assume that the tensor DWjk(p) is approximated over a given cell Om using Eq. (11) as follows: Z Z Nijk ðs; pÞDW jk ðpÞ dO ¼ N ijk ðs; pÞf‘ ðpÞ dO DW ‘jk (14)
3 X2 r
1 ξ1
θ S
X1 Fig. 3. Local coordinates.
The shape function Fk can also be expressed in terms of a cylindrical coordinate system with the origin at the collocation point s (Fig. 3).
Fk ðpÞ ¼ Fk ðsÞ þ
rðs; pÞ ½mk cosðyÞ þ nk senðyÞ 2A
(13)
b
(11)
Thus, the integral on the right hand side of Eq. (8) can be carried out over r and then transformed to a cell boundary integral as follows: Z Z 3 1 r 3 Fk ðpÞw ðQ ; pÞ dO ¼ Fk ðsÞ ln r 8 4 p D 4 Gm Om
r4 7 1 @r þ ½mk cos y þ nk sen y dG ln r 10 2A @n 5
(12) where Gm represents the cell contour formed by three side elements. Using this transformation one can evaluate the integral over the cell by performing three cell side integrals, as they are cell boundary elements. Therefore, similar to the boundary elements, the integrals over each cell side or cell boundary element can be computed either using the analytical expression or an appropriate numerical scheme to carry out the integrals over the boundary elements. Then, after performing the integral over all cell sides, the obtained coefficients are cumulated appropriately in a matrix Sbw (Eq. (13)). This matrix gives the contribution of DTNw that
Om
Om
where DWjk‘ represents the nodal values of the condensate incremental tensorDWjk‘(p). After carrying out the integral over the dimension r, Eq. (14) becomes: Z N ijk ðs; pÞDW jk ðpÞ dOðpÞ Om Z ff ðsÞ þ rðs; pÞ ½m‘ cosðyÞ þ n‘ senðyÞgr;n dG DW l (15) ¼ N ‘ ijk jk 4A Gm ¯ *ijk ¼ rN*ijk is given in the Appendix A. where r,n ¼ qr/qn and N Thus, one can perform the integral along the boundary element and along the cell sides as shown in Eq. (15) to obtain the algebraic form of Eq. (4) duly modified to be written in increments. Writing this algebraic equation for the appropriate number of collocation points, one obtains H us DU s ¼ Gus DP s þ S us DW yy þ DBus
(16)
where DUs and DPs contain all incremental boundary displacement and traction nodal values of the stretching problem, DWyy represents the increment of the rotation product defined at each domain node, Ssu is the corresponding matrix obtained by integrating all domain cells and DBsu gives the in-plane domain incremental load effects. Eqs. (13) and (16) are boundary algebraic relations obtained for collocation points defined along the boundary (or outside the
ARTICLE IN PRESS L. Waidemam, W.S. Venturini / Engineering Analysis with Boundary Elements 33 (2009) 1223–1230
q
q
q
q
x2
1227
x1
x2
x1
a
a
a
a
Fig. 4. Square plate with uniform load; (a) geometry; (b) discretization.
domain). To complete the system of equations for this problem, we have to obtain the algebraic relations from the integral equations of internal values. Thus, the integral Eqs. (5)–(7), transformed into algebraic relations, are the algebraic equations required to solve the problem. We have to write them in their rate forms, carry out the time integration to obtain the corresponding incremental representations, and then, by assuming boundary and domain discretizations and approximations already adopted to obtain Eqs. (13) and (16), they are transformed to
Dh ¼ H yb DU b þ Gyb DP b þ S yb DT wN þ DByb
(17)
DN ¼ ANs DU b þ DF Ns þ S Ns DW yy
(24)
y
where the matrix Ax is conveniently built using columns of Hxy and Gxy corresponding to unknown boundary values, and DFxy is a vector collecting all contributions of prescribed domain and boundary values. Solving Eqs. (20) and (21) leads to w w DX b ¼ DM w b þ Q b DT N
(25)
DX s ¼ DM us þ Q us DW yy y
[Axy]1
(26) y
Qxy
y 1
Sxy.
Dv ¼ Hwb DU b þ Gwb DP b þ S wb DT wN þ DBwb
(18)
with DMx ¼ DFx and ¼ [Ax ] Replacing Eqs. (25) and (26) accordingly into Eqs. (22)–(24) gives
DN ¼ H Ns DU s þ GNs DP s þ S Ns DW yy þ DBNs
(19)
Dh ¼ DN hb þ Q hb DT vN
(27)
Dv ¼ DN vb þ Q vb DT vN
(28)
where Dh, Dv and DN are vectors containing the rotation, curvature and membrane internal force increments at the domain nodes defined by the adopted discretization. Again, the domain integral terms in Eqs. (5)–(7) can be transformed to cell side integrals using a transformation similar to the one employed to obtain Eqs. (12) and (15). The integrals containing the kernels w*,i and w*,ij, in Eqs. (5) and (6), are not singular and lead to regular cell side integrals. The singular domain integral, in Eq. (7), leads to a 1/r singular integral along the cell sides. This integral is the same one already treated in nonlinear BEM formulations [29,32,33]. The sub-element scheme used to evaluate the element boundary element integrals has already been demonstrated to be accurate. It was verified by comparing the computed values with values obtained using analytical integrals. The properties of the resulting matrix SsN also demonstrated the accuracy of the adopted integration scheme. By applying a constant field DW we could verify that the global errors are less than 1010. After applying the boundary conditions, Eqs. (13) and (16)–(19) become: w b w Aw b DX ¼ DF b þ S b DT N
w
(20)
Aus DX s ¼ DF us þ S us DW yy
(21)
Dh ¼ Ayb DX b þ DF yb þ S yb DT wN
(22)
Dv ¼ Awb DX b þ DF wb þ S wb DT wN
(23)
DN ¼ DN Ns þ Q Ns DW yy y
y
(29) z
y
y
y
Ax Qxz+Sxy,
where DNx ¼ Ax DMx +DFx and Qx ¼ with z given by the superscript of the corresponding boundary algebraic equation used to compute Qxz, i.e. Qxw from Eq. (25) or Qxu from Eq. (26). Eqs. (25)–(29) are the necessary relations to solve a geometric nonlinear plate problem. However, one has to correctly treat the increments DWij(p) and DT(p). We can first find the rates of the densities in Eqs. (3) and (4) and obtain the incremental forms of DT and DW for a given time interval Dtn, as follows:
DT ¼ DN v þ N Dv
(30)
DW ¼ Dh h þ h Dh
(31)
4. System solution Eqs. (27)–(29) represent the nonlinear system to be solved. They represent the equilibrium at each internal point. The boundary conditions and the corresponding unknown boundary values are already taken into account. The increments DNxy are the elastic solution obtained after applying the load increment. The equation will be solved in terms of the increments Dh, Dv and DN. The increments DT and DW will be written in terms of the three chosen values according to Eqs. (30) and (31), as
ARTICLE IN PRESS 1228
L. Waidemam, W.S. Venturini / Engineering Analysis with Boundary Elements 33 (2009) 1223–1230
follows: y
y
y
F y ðDh; Dv; DNÞ ¼ Dh þ DN b þ Q b DN v þ Q b N Dv ¼ 0
(32)
F w ðDh; Dv; DNÞ ¼ Dv þ DN b þ Q b DN v þ Q b N Dv
w
w
w
(33)
N N F N ðDh; Dv; DNÞ ¼ DN þ DN N s þ Q s Dh h þ Q s h Dh
(34)
The above nonlinear system of equations is solved by applying the Newton–Raphson’s scheme. Within a time increment Dt ¼ tn+1tn an iterative process may be required to achieve the equilibrium. Then, from the solution at iteration i, the next attempt at iteration (i+1) may require additive increment corrections, i.e.
DN iþ1 n
¼
DN in
þ
dDN in
Remark: the C matrix above is sparse, and therefore Eq. (38) can be appropriately solved by using sparse system procedures, which are much more efficient in comparison with other techniques based on system sub-structuring techniques.
5. Numerical applications Three examples have been analysed using the proposed model to demonstrate the performance of the proposed formulation. Square plate subjected to a uniform load: the square plate chosen to run this example is defined by the side length equal to a, 2.5
(35) (36)
Dviþ1 ¼ Dvin þ dDvin n
(37)
(wmax /h)
2.0
Dhiþ1 ¼ Dhin þ dDhin n
To solve a nonlinear system of equations (see Simo and Hughes [34]), Taylor’s expansion is applied to the equilibrium Eqs. (32)–(34). We can assume that considering only the first term of this expansion is enough accurate to treat the proposed nonlinear system of equations. Thus, the equilibrium equations can be approached by 2 3 @F y @F y @F y i i i 7 9 8 9 6 6 @Dhn @Dvn @DN n 78 F ðDhin ; Dvin ; DN in Þ > > 6 7> dDhin > > > > > < = < y = 6 @F w 7 @F @F w w 7 dDvin þ ¼ 0 F w ðDhin ; Dvin ; DN in Þ þ 6 6 i i i 7 @ D v @ D N > > > 6 @Dhn n n 7> > : dDN i > ; : F ðDhi ; Dvi ; DN i Þ > ; 6 7> 6 @F N n N n n n @F N @F N 7 4 5 i i i @Dhn @Dvn @DN n
1.5 1.0 0.5 0.0 20
0
Linear solution
n
n
ME
IE
x2
a
a Fig. 7. Boundary and internal discretizations.
1.0 0.8 0.6 0.4 0.2 0.0 10
15
20
(q.a4/E.h4) Linear solution IE
Ye & Liu [14] - ME
q
(39)
5
Ye & Liu [14] - IE
x1
1.2
0
100
q
n
(wmax/h)
N
80
Fig. 6. Clamped plate—load displacement curve.
Thus, the corrections to be cumulated during the iterations are obtained by solving Eq. (38) 8 8 9 9 i F ðDhin ; Dvin ; DN in Þ > > > > > > > < y < dDhn > = = dDvin ¼ ½C1 F w ðDhin ; Dvin ; DN in Þ (40) > > > > > > : F ðDhi ; Dvi ; DN i Þ > : dDN i > ; ; n
60 (q.a4/16.D.h)
(38) where each term of the tangent matrix is given by 3 2 Q yb II win I Q yb N in II 7 6 w w 0 II þ Q b N in II Q b II vin 7 ½C ¼ 6 5 4 i i QN ðI h þ h IÞ 0 II s n n
40
Experimental [14]
Ye & Liu [14] - IE
ME
Ye & Liu [14] - ME
Fig. 5. Simply supported plate—load displacement curve.
25
120
ARTICLE IN PRESS L. Waidemam, W.S. Venturini / Engineering Analysis with Boundary Elements 33 (2009) 1223–1230
2.0 1.5 1.0 0.5 0.0 0
1
2
3
4
5
6
7
8
9
10
(q.a4/E.h) Linear solution
Present study
Wen et al [17]
Exact
Fig. 10. Load deflection curve for the circular plate.
Circular Clamped side: the third example is a circular plate also taken from the literature [17,15]. The geometric data and adopted discretization are given in Fig. 9. Forty linear elements were used along the whole boundary, while the domain was discretized using 100 internal cells. Poisson’s ratio is equal to 0.3, while the ratio h/a ¼ 0.02. The results obtained in terms of load deflection are given in Fig. 10. Good agreement is observed when compared with the analytical solution and the numerical results [17].
6. Conclusions A boundary element formulation to analyse von Ka´rma´n plates is proposed herein. The domain approximations were simplified by reducing the density of domain integrals into simple equivalent values. This scheme significantly reduces the computational effort needed to compute the matrices related to the domain values. The tangent operator has been derived, and the Newton process has been implemented to solve the nonlinear set of equations, which has lead to accurate and reliable computed values using a reduced number of iterations.
1.2 1.0 0.8 (wmax/h)
2.5
(wmax/h)
thickness h, the ratio h/a ¼ 0.01 and Poisson’s ratio v ¼ 0.3. The uniform applied load is q (Fig. 4). Several boundary conditions have been analysed. We started by assuming simple supported and clamped plate conditions. In each case, the in-plane boundary displacements can also be set equal to zero (IE) or kept free (ME). The obtained results, compared with other numerical and analytical solutions, demonstrate that the proposed formulation is accurate. Running several other internal and boundary meshes has also demonstrated that the convergence is obtained quickly. Even with a rather coarse mesh, particularly to integrate the domain integrals, the results are already accurate. The presented results were obtained by using 160 boundary elements and only 8 domain cells. The results obtained are shown in Figs. 5 and 6, for the simply supported plate and clamped plate respectively. As can be seen, the results are in total agreement with those given by Ye and Liu [14]. Square plate with two clamped sides; the geometry of this example is taken from the first application, and only the boundary conditions are changed. For the two opposite clamped sides, indicated in Fig. 6, the following values are prescribed: w ¼ 0, qw/qn ¼ 0, ux ¼ 0 and uy ¼ 0, while for the remaining simply supported sides the prescribed values are: w ¼ 0, Mn ¼ 0, ux ¼ 0 and uy ¼ 0. The boundary and internal discretization are now given by 80 boundary elements and 32 cells (Fig. 7). The ratio between the thickness and the side length was set to h/a ¼ 0.01. The computed load deflection curve obtained for this example is given in Fig. 8. As can be seen, the computed results are almost the same as those obtained by Purbolaksono and Aliabadi [18] using a BEM bending plate model based on Reissner’s theory.
1229
0.6 0.4
Acknowledgment
0.2 0.0 0
10
Linear solution
20
30 (q.a4/E.h4) Present study
40
50
To FAPESP—Sa˜o Paulo State Research Foundation for the support given to this work.
Purbolaksono & Aliabadi [18]
Appendix A Fig. 8. Square plate deflections with two opposite clamped sides.
Fundamental values used throughout the text: (a) For the plate bending problem: w ðs; pÞ ¼
1 2 r ðln r 1=2Þ 8pD
@w r ðs; pÞ ¼ ln rðr;k nk Þ @n 4pD
M ns ðs; pÞ ¼
2a Fig. 9. Geometry and the adopted discretization for the circular plate.
V n ðs; pÞ ¼
ð1 nÞ ðr;k sk Þðr;l nl Þ 4p
h i ð1 nÞ 1 ðr;k nk Þ 2ð1 nÞðr;l sl Þ2 3 þ n þ ðr;k nk Þ2 4pr 4pR
ARTICLE IN PRESS 1230
L. Waidemam, W.S. Venturini / Engineering Analysis with Boundary Elements 33 (2009) 1223–1230
(b) For the plate-stretching problem: uij ðs; pÞ ¼
1 ½ð3 4nÞdij lnðrÞ þ r;i r;j 8pGhð1 nÞ
pij ðs; pÞ ¼
1 @r ½ð1 2nÞdij þ 2r;i r;j ð1 2nÞðr;i nj r;j ni Þ 4pð1 nÞr @n
N ijk ðs; qÞ ¼
Dijk ðs; qÞ ¼
Sijk ðs; qÞ ¼
1 ½2r;i r;j r;k þ ð1 2n0 Þðdij r;k þ dik r;j djk r;i Þ 4pð1 n0 Þr
1 fð1 2nÞ½dik r;j þ djk r;i dij r;k þ 2r;i r;j r;k g 4prð1 nÞ
Gh 2p
r 2 ð1
nÞ
2
@r ð1 2nÞdij r;k þ nðdik r;j þ djk r;i Þ 4r;i r;j r;k @n
þ 2nðni r;j r;k Þ þ nj r;i r;k þ ð1 2nÞ ð2nk r;i r;j þ nj dik þ ni djk Þ ð1 4nÞnk dij g
with
n0 ¼ n=ð1 þ nÞ
References [1] Jaswon MA, Maiti M, Symm GT. Numerical biharmonic analysis and some applications. International Journal of Solids and Structures 1967;3:309–32. [2] Jaswon MA, Maiti M. An integral formulation of plate bending problems. Journal of Engineering Mathematics 1968;2:83–93. [3] Bezine GP. Boundary integral formulation for plate flexure with arbitrary boundary conditions. Mechanics Research Communications 1978;5(4): 197–206. [4] Tottenhan H. The boundary element method for plates and shells. In: Banerjee PK, Butterfield R., editors. Developments in boundary element methods, 1979, p. 173–205. [5] Stern MA. A general boundary integral formulation for the numerical solution of plate bending problems. International Journal of Solids and Structures 1979;15:769–82. [6] Van der Weee¨n F. Application of the boundary integral equation method to Reissner’s plate model. International Journal for Numerical Methods in Engineering 1982;18:1–10. [7] Barcellos CA, Silva LHM. A boundary element formulation for Mindlin’s plate model. In: Brebbia CA, Venturini WS, editors. BETECH-87. Southampton: CMP; 1897. p. 123–30. [8] Aliabadi MH. Plate bending analysis with boundary element. Southampton: Computational Mechanics Publications; 1998. [9] Moshaiov A, Vorus WS. Elastoplastic plate bending analysis by a boundary element method with initial plastic moments. International Journal of Solids and Structures 1986;22(11):1213–29. [10] Chueiri LHM, Venturini WS. Elastoplastic BEM to model concrete slabs. In: Brebbia CA, et al., editors. Boundary elements XVII. CMP; 1995. p. 149–56.
[11] Providakis CP. Quasi-static analysis of viscoplastic plates by the boundaryelement method. Engineering Analysis with Boundary Elements 1993;11(4): 265–8. [12] Morjaria M. Inelastic analysis of transverse deflection of plates by the boundary element method. Journal of Applied Mechanics—Transactions of the ASME 1980;47:291. [13] Kamiya N, Sawaki Y. An integral equation approach to finite deflection of elastic plates. International Journal of Non-Linear Mechanics 1982;17(3): 187–94. [14] Ye T-Q, Liu Y-J. Finite deflection analysis of elastic plate by the boundary element method. Applied Mathematical Modelling 1985;9:183–8. [15] Tanaka M, Matsumoto T, Zheng Z-D. Incremental analysis of finite deflection of elastic plates via boundary-domain-element method. Engineering Analysis with Boundary Elements 1996;17:123–31. [16] Wang W, Ji X, Tanaka M. A dual boundary element approach for the problems of large deflection of thin elastic plates. Computational Mechanics 2000; 26:58–65. [17] Wen PH, Aliabadi MH, Young A. Large deflection analysis of Reissner plate by boundary element method. Computers & Structures 2005;83(10–11):870–9. [18] Purbolaksono J, Aliabadi MH. Large deformation of shear-deformable plates by the boundary element method. Journal of Engineering Mathematics 2005;51(3):211–30. [19] Dirgantara T, Aliabadi MH. A boundary element formulation for geometrically nonlinear analysis of shear deformable shells. Computer Methods in Applied Mechanics and Engineering 2006;195(33–36):4635–54. [20] Supriyono, Aliabadi MH. Boundary element method for shear deformable plates with combined geometric and material nonlinearities. Engineering Analysis with Boundary Elements 2006;30(1):31–42. [21] Fung YC. Foundations of solid mechanics. Englewood Cliff, NJ: Prentice-Hall, Inc.; 1964. [22] Brebbia CA, Telles JCF, Wrobel LC. Boundary element techniques. Theory and applications in engineering. Berlin and New York: Springer-Verlag; 1984. [23] Bui HD. Some remarks about the formation of three-dimensional thermoelastoplastic problems by integral equations. International Journal of Solids and Structures 1978;14:935–9. [24] Venturini WS, Paiva JB. Boundary elements for plate bending analysis. Engineering Analysis with Boundary Elements 1993;11:1–8. [25] Paiva JB, Venturini WS. Alternative technique for the solution of plate bending problems using the boundary element method. Advances in Engineering Software 1992;14(4):265–71. [26] Fernandes GR, Chaves EWV, Venturini WS. Plate bending boundary element formulation considering variable thickness. Engineering Analysis with Boundary Elements 1999;23:405–18. [27] Leite LGS, Venturini WS. Stiff and soft thin inclusions in two-dimensional solids by the boundary element method. Engineering Analysis with Boundary Elements 2003;29(3):257–67. [28] Leite LGS, Venturini WS. Stiff and soft thin inclusions in two-dimensional solids by the boundary element method. Engineering Analysis with Boundary Elements 2005;29(3):257–67. [29] Botta AS, Venturini WS, Benallal A. BEM applied to damage models emphasizing localization and associated regularization techniques. Engineering Analysis with Boundary Elements 2005;29:814–27. [30] Fernandes GR, Venturini WS. Building floor analysis by the boundary element method. Computational Mechanics 2005;35:277–91. [31] Fernandes GR, Venturini WS. Non-linear boundary element analysis of floor slabs reinforced with rectangular beams. Engineering Analysis with Boundary Elements 2007;31:721–37. [32] Bonnet M, Mukherjee S. Implicit BEM formulation for usual and sensitivity problems in elasto-plasticity using the consistent tangent operator concept. International Journal of Solids and Structures 1996;33:4461–80. [33] Benallal A, Fudoli CA, Venturini WS. An implicit BEM formulation for Gradient plasticity and localization phenomena. International Journal for Numerical Methods in Engineering 2002;53:1853–69. [34] Simo JC, Hughes TJR. Computational inelasticity. New York: Springer; 1998.
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
BEM formulation for von Ka´rma´n plates Leandro Waidemam a, Wilson S. Venturini b, a b
´ Federal University of Technology, BR 369, km 0.5, 87301-006-Campo Moura ˜o, Brazil Parana ˜o Carlos School of Engineering, University of Sa ˜o Paulo, Av. Trabalhador Sa ˜o-Carlense 400, 13566-970 Sa ˜o Carlos, Brazil Sa
a r t i c l e in fo
abstract
Article history: Received 1 April 2007 Accepted 24 April 2009 Available online 7 June 2009
This work deals with nonlinear geometric plates in the context of von Ka´rma´n’s theory. The formulation is written such that only the boundary in-plane displacement and deflection integral equations for boundary collocations are required. At internal points, only out-of-plane rotation, curvature and in-plane internal force representations are used. Thus, only integral representations of these values are derived. The nonlinear system of equations is derived by approximating all densities in the domain integrals as single values, which therefore reduces the computational effort needed to evaluate the domain value influences. Hyper-singular equations are avoided by approximating the domain values using only internal nodes. The solution is obtained using a Newton scheme for which a consistent tangent operator was derived. & 2009 Elsevier Ltd. All rights reserved.
Keywords: Bending plates Geometrical nonlinearities
1. Introduction The boundary element method (BEM) has been successfully applied to solve plate bending problems many times. An important characteristic of the boundary methods applied to plate bending is that all boundary values are approximated by the same shape function, which avoids the necessity of using higher order derivatives of the displacement approximation to compute internal forces. Thus, bending and twisting moments and shear forces are precisely evaluated. The method has already proved to be accurate and reliable enough for this kind of application. The first works using the indirect version of the boundary methods applied to Kirchhoff’s plates were provided by Jaswon et al. [1] and Jaswon and Maiti [2]. The direct formulation based on the same simplified plate bending hypothesis was proposed by Bezine [3], Tottenhan [4] and Stern [5]. Van der Weee¨n [6] has derived the boundary element formulation in the context of the Reissner’s theory of plates. Barcellos and Silva proposed the extension to consider the Mindlin hypothesis [7]. More recent contributions in boundary elements to treat shear deformable plates can be found in Aliabadi [8]. The plate bending numerical formulation is a very important subject in engineering because of its many applications to a large number of complex structures, such as aircraft, ships, cars, pressure vessels, and off shore structures, among others. Usually, these complex problems require accurate plate bending models to
Corresponding author.
E-mail addresses: [email protected] (L. Waidemam), [email protected] (W.S. Venturini). 0955-7997/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2009.04.012
take into account physical and geometric nonlinear effects. Regarding the use of plastic models in conjunction with plate bending BEM models, we have to point out the works carried out by Moshaiov and Vorus [9]. Chueiri and Venturini [10] who have proposed a general model that takes into account not only the plastic criterion written in terms of moments, by also the integration of the corrector stress field over the thickness, which is very convenient for the modelling of reinforced concrete behaviour. Providakis [11] was the first to propose a viscoplastic model for plate bending analysis. Several BEM formulations have already been proposed to consider nonlinear geometric effects in the global equilibrium of plate bending problems. One of the first works treating this subject was proposed by Morjaria [12]. Kamiya and Sawaki [13] have proposed a BEM formulation for elastic plates governed by the Berger equation. The first BEM formulation to analyse plate bending problems within the context of von Ka´rma´n hypothesis was proposed by Ye and Liu [14], who used a fictitious load distributed over the domain to model nonlinear effects. The von Ka´rma´n hypothesis was also adopted by Tanaka et al. [15] to develop a more elaborated BEM incremental formulation to deal with finite deflections of thin elastic plates. Wang et al. [16] have worked on von Ka´rma´n plates and introduced the dual reciprocity approach based on global radial functions to approximate the correcting integral term. All works reported above were proposed in the context of thin plates. Several other important works have appeared more recently that demonstrate the efficiency with which BEM formulations deal with shear deformable plates based on the Reissner–Mindlin hypothesis and taking into account geometric and material nonlinearities. Among them, we emphasise the works carried out by Aliabadi and his
ARTICLE IN PRESS 1224
L. Waidemam, W.S. Venturini / Engineering Analysis with Boundary Elements 33 (2009) 1223–1230
co-workers: Wen, et al. [17], Purbolaksono and Aliabadi [18], Dirgantara and Aliabadi [19], and Supriyono and Aliabadi [20]. In this work, we return to the BEM formulation based on von Ka´rma´n’s theory. Emphasis is placed on the accurate evaluation of the domain integrals that are approximated with cells and to the solution technique for which a tangent consistent operator is proposed. Examples of plates with finite deflection are analysed, and the results are compared with other numerical solutions.
þ
Nc X
Rc ðC k Þwc ðS; C k Þ þ
k¼1
Z þ
O
Z
gðpÞw ðS; pÞdO Og
N ij ðpÞw;ij ðpÞw ðS; pÞ dO Z
(3)
P ij ðS; PÞuj ðPÞ dG þ
Z
uij ðS; PÞpj ðPÞ dG Z Z 1 þ uij ðS; pÞbj ðpÞ dO N ðS; pÞw;j ðpÞw;k ðpÞ dO 2 O ijk O
C ij ðSÞuj ðSÞ ¼
G
G
(4)
2. Theoretical background Without a loss of generality, let us consider a single thin plate region O with boundary G over which a distributed load q is applied orthogonally to the middle surface, i.e. in the direction x3, as shown in Fig. 1. This plate region can also be subjected to in-plane forces (directions x1 and x2) that are either distributed over the domain or applied along the boundary. To write the field equations for this plate problem, the following hypothesis can be assumed according to the von Ka´rma´n theory: the strains are assumed to be small and the final deflections are on the order of the plate thickness h. For this problem, the differential equations of the bending and in-plane problems are given by [21] mij ;ij þ N ij w;ij bi w;i þ mi ;i þ q ¼ 0
(1a)
Nij ;j þ bi ¼ 0
(1b)
where w represents the deflections, mij are bending and twisting moments, Nij gives the in-plane internal forces, bi represents the in-plane domain loads applied and mi is the applied moment over the plate domain; the subscripts are in the range i,j ¼ {1, 2}. The internal forces in Eqs. (1a) and (1b) can be written in terms of displacements as follows: mij ¼ Dðndij w;kk þ ð1 nÞw;ij Þ
where ( )* represents the well-known fundamental solutions of both problems; Vn, Mn, w and w,n are the effective shear forces, moments, deflections and rotations orthogonal to the plate boundary applied along the boundary, respectively; Rc and wc represent the corner reactions and deflections, respectively; pj and uj represent, respectively, the in-plane tractions and displacements along the plate boundary; w,j, w,ij and Nij are the rotations, curvatures and in-plane forces, respectively; S, P, p and Ck represent the collocation point and the field points defined along the boundary, inside the domains and at the corners, respectively; C(S) and Cij(S) are the well-known free terms of the plate bending and stretching problems [22]. To complete the required integral equations, one has to obtain the integral representations of rotations, curvatures and internal forces for a given domain collocation s as follows: Z w;i ðsÞ ¼ fV n ;i ðs; PÞwðPÞ M n ;i ðs; PÞw;n ðPÞg dG G
Nc X k¼1
þ
Nc X
Eh ð1 nÞ
2
nðuk ;k þ w;k w;k =2Þdij þ
Z
(2a) 1n ðui ;j þ uj ;i þ w;j w;i Þ 2
w;ij ðsÞ ¼ (2b)
CðSÞwðSÞ þ
Z G
Z ¼ G
ðV n ðS; PÞwðPÞ M n ðS; PÞw;n ðPÞÞ dG þ
Nc X
Rc ðS; C k Þwc ðC k Þ
k¼1
ðV n ðPÞw ðS; PÞ Mn ðPÞw;n ðS; PÞÞ dG
Γ x1 x2
Ω
x3
Fig. 1. General plate domain.
Z G
Nc X
Rc ;ij ðs; C k Þwc ðC k Þ þ
Nc X
Rc ðC k Þwc ;ij ðs; C k Þ þ
k¼1
Z þ
O
Nij ðsÞ ¼
G
fV n ðPÞw;i ðs; PÞ M n ðPÞw;ni ðs; PÞg dG
Z Og
gðpÞw;i ðs; pÞ dO
(5)
fV n ;ij ðs; PÞwðPÞ M n ;ij ðs; PÞw;n ðPÞg dG
k¼1
þ
Z
N jk ðpÞw;jk ðpÞw;i ðs; pÞ dO
þ
where ui are the in-plane displacement components. The problem definition is then completed by assuming the following boundary conditions over G: ui ¼ u% i on G1 (generalised displacements, in-plane displacements, deflections and rotations); and pi ¼ p¯i on G2 (generalised tractions, in-plane tractions normal bending moment and effective shear forces), where G1[G2 ¼ G. To obtain the integral equations of the plate bending and stretching problems considering geometrical nonlinearities within the context previously defined, one can apply the Betti’s reciprocity to the linear parts of the stress and strain fields. The integral representations of both problems are given by [22]:
Rc ðC k Þwc ;i ðs; C k Þ þ
k¼1
O
Nij ¼
Rc ;i ðs; C k Þwc ðC k Þ þ
Z G
fV n ðPÞw;ij ðs; PÞ M n ðPÞw;nij ðs; PÞg dG
Z Og
gðpÞw;ij ðs; pÞ dO
N km ðpÞw;km ðpÞw;ij ðs; pÞ dO
Z Dijk ðs; PÞpk ðpÞ dG Sijk ðs; PÞuk ðPÞ dG G G Z Z 1 þ Dijk ðs; pÞbk ðpÞ dO T ðs; pÞw;k ðpÞw;l ðpÞ dO 2 O ijkl O Gh ð2w;i ðsÞw;j ðsÞ þ w;k ðsÞw;k ðsÞdij Þ þ 8ð1 u Þ
(6)
Z
(7)
where v¯ ¼ v/(1+v) is used to simulate the plane stress conditions; the higher order tensors are given in the appendix. The free term appearing in Eq. (7) is due to the differentiation of a singular integral. It has been derived for 3D, plane stress and plane strain problems by Bui [23] and extended for several other integral equations containing kernels with the same singularity.
Ωg 3. Algebraic equations Eqs. (3)–(7) were obtained from von Ka´rma´n’s field equations; therefore, they are written in terms of total values. To reduce the number of domain integrals in this section, it is convenient to replace the densities given by the product of the two values,
ARTICLE IN PRESS L. Waidemam, W.S. Venturini / Engineering Analysis with Boundary Elements 33 (2009) 1223–1230
which appear in the domain integrals, by single equivalent values. Then, the new equivalent values will be approximated by using node values and the shaped functions. Thus, the domain value Nij(p)w,ij(p) at a field point p that appears in the bending problem Eqs. (3), (5) and (6) will be replaced by an equivalent scalar value T(p). Similarly, in Eqs. (4) and (7) of the stretching problem, we replace the domain density w,j(p)w,k(p) by an equivalent tensor Wjk(p). As is typical when solving a nonlinear problem using an incremental process, it is convenient to write the integral representations derived in the previous section into their corresponding rate integral representations. Bearing in mind that we have already replaced the densities in the bending problem equations by the equivalent scalar T(p) and in the in-plane problem equations by the equivalent tensor Wjk(p), we now have to ˙ jk(p), replace these values by the corresponding rates T˙(p) and W respectively. All the other boundary and domain values appearing in all integral equations derived in the previous sections are also replaced by their rates. The starting point of the numerical formulation for the formulated boundary value problem is its incremental representation. The governing equations and the derived integral representations have to be written in increments. To find their appropriate incremental forms, we assume that Dt ¼ tn+1tn is a typical time-step in the time discretization. Then, any rate quantity x˙ integrated along this time interval Dt becomes Dx ¼ xn+1xn, which will replace x˙ in all governing equations derived from the integral representations. Thus, the rate integral representations (3)–(7) are transformed into incremental representations. For instance, Eq. (3) will be first written in rates and then in increments as follows: CðSÞDwðSÞ þ
Z G
þ
Nc X
ðV n ðS; PÞDwðPÞ Mn ðS; PÞDw;n ðPÞÞ dG
Rc ðS; C k ÞDwc ðC k Þ
k¼1
Z ¼
G
þ
ðDV n ðPÞw ðS; PÞ DM n ðPÞw;n ðS; PÞÞ dG
Nc X
DRc ðC k Þwc ðS; C k Þ þ
k¼1
Z
DTðpÞw ðS; pÞ dO
þ
Z
DgðpÞw ðS; pÞ dO
Og
(8)
O
The same modification has to be performed on the other integral representations, Eqs. (4)–(7). After obtaining the finite step integral equations, one has to transform them into algebraic equations by approximating the boundary and domain incremental values. For this purpose, the plate boundary G is discretized into a series of elements, Gs, along which generalised displacements and tractions are approximated by using continuous and discontinuous linear boundary elements. For the present plate stretching–bending problem, discontinuities of boundary values are always expected, particularly at corners and traction jumps. The boundary value discontinuities are always introduced by defining the collocation points placed inside the elements. Typically, for the plate stretching and bending boundary element formulations, we need two equations for each problem for each boundary node used to discretize the boundary. Thus, two algebraic equations for each node are referred to the stretching problems and derived from the incremental form of Eq. (4). Another two algebraic equations for the plate bending problem are derived from the deflection Eq. (8). Since we are avoiding writing the slope equation for this purpose, we have to write two independent deflection integral representations for
1225
each node; therefore, two collocations for each boundary node have to be defined: one on the boundary, either defined at the node or inside the element, and the second out of the domain. Moreover, the corner reactions can also be properly accounted for in the effective shear forces of the adjacent boundary elements, avoiding an additional unknown corner value. This scheme has proved to lead to stable and accurate results. As previously demonstrated, the required outside collocations have to be defined at points very near the boundary to guarantee the accuracy and stability of the solution [24,25]. For this work, they have been defined within the range of 0 to half of the adjacent element average length. It is important to stress that corner double nodes are always used to guarantee the discontinuities of all boundary values. Thus, the boundary values that refer to nodes defined at both elements adjacent to a corner are always independent of each other. Moreover, extra values of the deflection and the concentrated reaction can be associated with the corner. If these values are assumed to be independent, one extra equation for this corner has to be written. In this work, we prefer to replace the corner deflection by the average deflection values of the two double nodes, while the concentrated reaction is properly distributed along the adjacent elements and added to the effective shear force [26]. Accurately performing the integrals over the elements requires special care, particularly when carrying out quasi-singular integrals. To assure the accuracy of these integrals, an appropriate sub-element scheme was employed. This sub-element scheme has already been successfully used to compute quasi-singular integrals that result from the displacement and stress equations for linear and nonlinear stretching problems [27–29]. The accuracy of this scheme has also been confirmed when applied to plate bending problems [30,31]. In these works, the accuracy of the adopted scheme that performs numerically quasi-singular integrals was certified by comparing the results with the values obtained by using analytical expressions and also by verifying the rigid body movements. To obtain the algebraic forms of the Eqs. (4)–(8) the domain is also discretized into continuous and discontinuous linear cells over which the relevant domain values are approximated, wherever required. Therefore, the domain integrals are all approximated over the cells using the appropriate shape functions and nodal values. For these approximations, only domain values are used to approach the domain fields; therefore, discontinuous cells are always required for cells adjacent to the boundary (Fig. 2). Moreover, boundary and domain discretizations may be completely independent of each other. Defining cells with sides larger than the corresponding adjacent boundary elements is recommended to avoid choosing internal nodes too close to the boundary. To describe the scheme used to compute the integrals over a given cell Om, let us take the last domain integral in Eq. (8). Assuming the density DT(p) is approximated over the given cell Om using the shape function Fk(p) and the nodal values DTwk, one can write Z
DTðpÞw ðs; PÞ dO ¼
Om
Z Om
Fk ðpÞw ðs; PÞ dODT kw
(9)
For the linear approximation case, the shape function Fk at a given field point p, written in terms of a global Cartesian system of coordinates, is
Fm ðpÞ ¼
1 ½am þ mm xðpÞ þ nm yðpÞ 2A
(10)
where A is the cell area and the constants am, mm and nm are the well-known coefficients of the linear approximation [22].
ARTICLE IN PRESS 1226
L. Waidemam, W.S. Venturini / Engineering Analysis with Boundary Elements 33 (2009) 1223–1230
Fig. 2. Linear domain approximation: (a) over continuous cells; (b) over discontinuous cells.
represents the increment of the product of the in-plane forces by curvatures. Although the integrals along the cell sides in Eq. (12) are not singular, one can use the same proposed sub-element scheme to compute the quasi-singular boundary elements to improve the computed vales. After performing all boundary and cell integration and storing the coefficients into the corresponding matrices, Eq. (8) becomes
Γ ξ2
Ω
2
w
w w b b w Hw b DU ¼ Gb DP þ S b DT N þ DBb
p
b
where DU and DP contain all boundary displacement and traction nodal values of the bending problem, plus the additional corner displacements and reactions, DTNw represents the condensed summation of the in-plane forces multiplied by the curvatures, Sbw is the corresponding matrix obtained by integrating all domain cells and DBbw gives the domain load effects. Similarly, one can also transform the in-plane displacement integral representation, Eq. (4), which has already modified by replacing the boundary and internal values by their increments, into its algebraic form. We first assume that the tensor DWjk(p) is approximated over a given cell Om using Eq. (11) as follows: Z Z Nijk ðs; pÞDW jk ðpÞ dO ¼ N ijk ðs; pÞf‘ ðpÞ dO DW ‘jk (14)
3 X2 r
1 ξ1
θ S
X1 Fig. 3. Local coordinates.
The shape function Fk can also be expressed in terms of a cylindrical coordinate system with the origin at the collocation point s (Fig. 3).
Fk ðpÞ ¼ Fk ðsÞ þ
rðs; pÞ ½mk cosðyÞ þ nk senðyÞ 2A
(13)
b
(11)
Thus, the integral on the right hand side of Eq. (8) can be carried out over r and then transformed to a cell boundary integral as follows: Z Z 3 1 r 3 Fk ðpÞw ðQ ; pÞ dO ¼ Fk ðsÞ ln r 8 4 p D 4 Gm Om
r4 7 1 @r þ ½mk cos y þ nk sen y dG ln r 10 2A @n 5
(12) where Gm represents the cell contour formed by three side elements. Using this transformation one can evaluate the integral over the cell by performing three cell side integrals, as they are cell boundary elements. Therefore, similar to the boundary elements, the integrals over each cell side or cell boundary element can be computed either using the analytical expression or an appropriate numerical scheme to carry out the integrals over the boundary elements. Then, after performing the integral over all cell sides, the obtained coefficients are cumulated appropriately in a matrix Sbw (Eq. (13)). This matrix gives the contribution of DTNw that
Om
Om
where DWjk‘ represents the nodal values of the condensate incremental tensorDWjk‘(p). After carrying out the integral over the dimension r, Eq. (14) becomes: Z N ijk ðs; pÞDW jk ðpÞ dOðpÞ Om Z ff ðsÞ þ rðs; pÞ ½m‘ cosðyÞ þ n‘ senðyÞgr;n dG DW l (15) ¼ N ‘ ijk jk 4A Gm ¯ *ijk ¼ rN*ijk is given in the Appendix A. where r,n ¼ qr/qn and N Thus, one can perform the integral along the boundary element and along the cell sides as shown in Eq. (15) to obtain the algebraic form of Eq. (4) duly modified to be written in increments. Writing this algebraic equation for the appropriate number of collocation points, one obtains H us DU s ¼ Gus DP s þ S us DW yy þ DBus
(16)
where DUs and DPs contain all incremental boundary displacement and traction nodal values of the stretching problem, DWyy represents the increment of the rotation product defined at each domain node, Ssu is the corresponding matrix obtained by integrating all domain cells and DBsu gives the in-plane domain incremental load effects. Eqs. (13) and (16) are boundary algebraic relations obtained for collocation points defined along the boundary (or outside the
ARTICLE IN PRESS L. Waidemam, W.S. Venturini / Engineering Analysis with Boundary Elements 33 (2009) 1223–1230
q
q
q
q
x2
1227
x1
x2
x1
a
a
a
a
Fig. 4. Square plate with uniform load; (a) geometry; (b) discretization.
domain). To complete the system of equations for this problem, we have to obtain the algebraic relations from the integral equations of internal values. Thus, the integral Eqs. (5)–(7), transformed into algebraic relations, are the algebraic equations required to solve the problem. We have to write them in their rate forms, carry out the time integration to obtain the corresponding incremental representations, and then, by assuming boundary and domain discretizations and approximations already adopted to obtain Eqs. (13) and (16), they are transformed to
Dh ¼ H yb DU b þ Gyb DP b þ S yb DT wN þ DByb
(17)
DN ¼ ANs DU b þ DF Ns þ S Ns DW yy
(24)
y
where the matrix Ax is conveniently built using columns of Hxy and Gxy corresponding to unknown boundary values, and DFxy is a vector collecting all contributions of prescribed domain and boundary values. Solving Eqs. (20) and (21) leads to w w DX b ¼ DM w b þ Q b DT N
(25)
DX s ¼ DM us þ Q us DW yy y
[Axy]1
(26) y
Qxy
y 1
Sxy.
Dv ¼ Hwb DU b þ Gwb DP b þ S wb DT wN þ DBwb
(18)
with DMx ¼ DFx and ¼ [Ax ] Replacing Eqs. (25) and (26) accordingly into Eqs. (22)–(24) gives
DN ¼ H Ns DU s þ GNs DP s þ S Ns DW yy þ DBNs
(19)
Dh ¼ DN hb þ Q hb DT vN
(27)
Dv ¼ DN vb þ Q vb DT vN
(28)
where Dh, Dv and DN are vectors containing the rotation, curvature and membrane internal force increments at the domain nodes defined by the adopted discretization. Again, the domain integral terms in Eqs. (5)–(7) can be transformed to cell side integrals using a transformation similar to the one employed to obtain Eqs. (12) and (15). The integrals containing the kernels w*,i and w*,ij, in Eqs. (5) and (6), are not singular and lead to regular cell side integrals. The singular domain integral, in Eq. (7), leads to a 1/r singular integral along the cell sides. This integral is the same one already treated in nonlinear BEM formulations [29,32,33]. The sub-element scheme used to evaluate the element boundary element integrals has already been demonstrated to be accurate. It was verified by comparing the computed values with values obtained using analytical integrals. The properties of the resulting matrix SsN also demonstrated the accuracy of the adopted integration scheme. By applying a constant field DW we could verify that the global errors are less than 1010. After applying the boundary conditions, Eqs. (13) and (16)–(19) become: w b w Aw b DX ¼ DF b þ S b DT N
w
(20)
Aus DX s ¼ DF us þ S us DW yy
(21)
Dh ¼ Ayb DX b þ DF yb þ S yb DT wN
(22)
Dv ¼ Awb DX b þ DF wb þ S wb DT wN
(23)
DN ¼ DN Ns þ Q Ns DW yy y
y
(29) z
y
y
y
Ax Qxz+Sxy,
where DNx ¼ Ax DMx +DFx and Qx ¼ with z given by the superscript of the corresponding boundary algebraic equation used to compute Qxz, i.e. Qxw from Eq. (25) or Qxu from Eq. (26). Eqs. (25)–(29) are the necessary relations to solve a geometric nonlinear plate problem. However, one has to correctly treat the increments DWij(p) and DT(p). We can first find the rates of the densities in Eqs. (3) and (4) and obtain the incremental forms of DT and DW for a given time interval Dtn, as follows:
DT ¼ DN v þ N Dv
(30)
DW ¼ Dh h þ h Dh
(31)
4. System solution Eqs. (27)–(29) represent the nonlinear system to be solved. They represent the equilibrium at each internal point. The boundary conditions and the corresponding unknown boundary values are already taken into account. The increments DNxy are the elastic solution obtained after applying the load increment. The equation will be solved in terms of the increments Dh, Dv and DN. The increments DT and DW will be written in terms of the three chosen values according to Eqs. (30) and (31), as
ARTICLE IN PRESS 1228
L. Waidemam, W.S. Venturini / Engineering Analysis with Boundary Elements 33 (2009) 1223–1230
follows: y
y
y
F y ðDh; Dv; DNÞ ¼ Dh þ DN b þ Q b DN v þ Q b N Dv ¼ 0
(32)
F w ðDh; Dv; DNÞ ¼ Dv þ DN b þ Q b DN v þ Q b N Dv
w
w
w
(33)
N N F N ðDh; Dv; DNÞ ¼ DN þ DN N s þ Q s Dh h þ Q s h Dh
(34)
The above nonlinear system of equations is solved by applying the Newton–Raphson’s scheme. Within a time increment Dt ¼ tn+1tn an iterative process may be required to achieve the equilibrium. Then, from the solution at iteration i, the next attempt at iteration (i+1) may require additive increment corrections, i.e.
DN iþ1 n
¼
DN in
þ
dDN in
Remark: the C matrix above is sparse, and therefore Eq. (38) can be appropriately solved by using sparse system procedures, which are much more efficient in comparison with other techniques based on system sub-structuring techniques.
5. Numerical applications Three examples have been analysed using the proposed model to demonstrate the performance of the proposed formulation. Square plate subjected to a uniform load: the square plate chosen to run this example is defined by the side length equal to a, 2.5
(35) (36)
Dviþ1 ¼ Dvin þ dDvin n
(37)
(wmax /h)
2.0
Dhiþ1 ¼ Dhin þ dDhin n
To solve a nonlinear system of equations (see Simo and Hughes [34]), Taylor’s expansion is applied to the equilibrium Eqs. (32)–(34). We can assume that considering only the first term of this expansion is enough accurate to treat the proposed nonlinear system of equations. Thus, the equilibrium equations can be approached by 2 3 @F y @F y @F y i i i 7 9 8 9 6 6 @Dhn @Dvn @DN n 78 F ðDhin ; Dvin ; DN in Þ > > 6 7> dDhin > > > > > < = < y = 6 @F w 7 @F @F w w 7 dDvin þ ¼ 0 F w ðDhin ; Dvin ; DN in Þ þ 6 6 i i i 7 @ D v @ D N > > > 6 @Dhn n n 7> > : dDN i > ; : F ðDhi ; Dvi ; DN i Þ > ; 6 7> 6 @F N n N n n n @F N @F N 7 4 5 i i i @Dhn @Dvn @DN n
1.5 1.0 0.5 0.0 20
0
Linear solution
n
n
ME
IE
x2
a
a Fig. 7. Boundary and internal discretizations.
1.0 0.8 0.6 0.4 0.2 0.0 10
15
20
(q.a4/E.h4) Linear solution IE
Ye & Liu [14] - ME
q
(39)
5
Ye & Liu [14] - IE
x1
1.2
0
100
q
n
(wmax/h)
N
80
Fig. 6. Clamped plate—load displacement curve.
Thus, the corrections to be cumulated during the iterations are obtained by solving Eq. (38) 8 8 9 9 i F ðDhin ; Dvin ; DN in Þ > > > > > > > < y < dDhn > = = dDvin ¼ ½C1 F w ðDhin ; Dvin ; DN in Þ (40) > > > > > > : F ðDhi ; Dvi ; DN i Þ > : dDN i > ; ; n
60 (q.a4/16.D.h)
(38) where each term of the tangent matrix is given by 3 2 Q yb II win I Q yb N in II 7 6 w w 0 II þ Q b N in II Q b II vin 7 ½C ¼ 6 5 4 i i QN ðI h þ h IÞ 0 II s n n
40
Experimental [14]
Ye & Liu [14] - IE
ME
Ye & Liu [14] - ME
Fig. 5. Simply supported plate—load displacement curve.
25
120
ARTICLE IN PRESS L. Waidemam, W.S. Venturini / Engineering Analysis with Boundary Elements 33 (2009) 1223–1230
2.0 1.5 1.0 0.5 0.0 0
1
2
3
4
5
6
7
8
9
10
(q.a4/E.h) Linear solution
Present study
Wen et al [17]
Exact
Fig. 10. Load deflection curve for the circular plate.
Circular Clamped side: the third example is a circular plate also taken from the literature [17,15]. The geometric data and adopted discretization are given in Fig. 9. Forty linear elements were used along the whole boundary, while the domain was discretized using 100 internal cells. Poisson’s ratio is equal to 0.3, while the ratio h/a ¼ 0.02. The results obtained in terms of load deflection are given in Fig. 10. Good agreement is observed when compared with the analytical solution and the numerical results [17].
6. Conclusions A boundary element formulation to analyse von Ka´rma´n plates is proposed herein. The domain approximations were simplified by reducing the density of domain integrals into simple equivalent values. This scheme significantly reduces the computational effort needed to compute the matrices related to the domain values. The tangent operator has been derived, and the Newton process has been implemented to solve the nonlinear set of equations, which has lead to accurate and reliable computed values using a reduced number of iterations.
1.2 1.0 0.8 (wmax/h)
2.5
(wmax/h)
thickness h, the ratio h/a ¼ 0.01 and Poisson’s ratio v ¼ 0.3. The uniform applied load is q (Fig. 4). Several boundary conditions have been analysed. We started by assuming simple supported and clamped plate conditions. In each case, the in-plane boundary displacements can also be set equal to zero (IE) or kept free (ME). The obtained results, compared with other numerical and analytical solutions, demonstrate that the proposed formulation is accurate. Running several other internal and boundary meshes has also demonstrated that the convergence is obtained quickly. Even with a rather coarse mesh, particularly to integrate the domain integrals, the results are already accurate. The presented results were obtained by using 160 boundary elements and only 8 domain cells. The results obtained are shown in Figs. 5 and 6, for the simply supported plate and clamped plate respectively. As can be seen, the results are in total agreement with those given by Ye and Liu [14]. Square plate with two clamped sides; the geometry of this example is taken from the first application, and only the boundary conditions are changed. For the two opposite clamped sides, indicated in Fig. 6, the following values are prescribed: w ¼ 0, qw/qn ¼ 0, ux ¼ 0 and uy ¼ 0, while for the remaining simply supported sides the prescribed values are: w ¼ 0, Mn ¼ 0, ux ¼ 0 and uy ¼ 0. The boundary and internal discretization are now given by 80 boundary elements and 32 cells (Fig. 7). The ratio between the thickness and the side length was set to h/a ¼ 0.01. The computed load deflection curve obtained for this example is given in Fig. 8. As can be seen, the computed results are almost the same as those obtained by Purbolaksono and Aliabadi [18] using a BEM bending plate model based on Reissner’s theory.
1229
0.6 0.4
Acknowledgment
0.2 0.0 0
10
Linear solution
20
30 (q.a4/E.h4) Present study
40
50
To FAPESP—Sa˜o Paulo State Research Foundation for the support given to this work.
Purbolaksono & Aliabadi [18]
Appendix A Fig. 8. Square plate deflections with two opposite clamped sides.
Fundamental values used throughout the text: (a) For the plate bending problem: w ðs; pÞ ¼
1 2 r ðln r 1=2Þ 8pD
@w r ðs; pÞ ¼ ln rðr;k nk Þ @n 4pD
M ns ðs; pÞ ¼
2a Fig. 9. Geometry and the adopted discretization for the circular plate.
V n ðs; pÞ ¼
ð1 nÞ ðr;k sk Þðr;l nl Þ 4p
h i ð1 nÞ 1 ðr;k nk Þ 2ð1 nÞðr;l sl Þ2 3 þ n þ ðr;k nk Þ2 4pr 4pR
ARTICLE IN PRESS 1230
L. Waidemam, W.S. Venturini / Engineering Analysis with Boundary Elements 33 (2009) 1223–1230
(b) For the plate-stretching problem: uij ðs; pÞ ¼
1 ½ð3 4nÞdij lnðrÞ þ r;i r;j 8pGhð1 nÞ
pij ðs; pÞ ¼
1 @r ½ð1 2nÞdij þ 2r;i r;j ð1 2nÞðr;i nj r;j ni Þ 4pð1 nÞr @n
N ijk ðs; qÞ ¼
Dijk ðs; qÞ ¼
Sijk ðs; qÞ ¼
1 ½2r;i r;j r;k þ ð1 2n0 Þðdij r;k þ dik r;j djk r;i Þ 4pð1 n0 Þr
1 fð1 2nÞ½dik r;j þ djk r;i dij r;k þ 2r;i r;j r;k g 4prð1 nÞ
Gh 2p
r 2 ð1
nÞ
2
@r ð1 2nÞdij r;k þ nðdik r;j þ djk r;i Þ 4r;i r;j r;k @n
þ 2nðni r;j r;k Þ þ nj r;i r;k þ ð1 2nÞ ð2nk r;i r;j þ nj dik þ ni djk Þ ð1 4nÞnk dij g
with
n0 ¼ n=ð1 þ nÞ
References [1] Jaswon MA, Maiti M, Symm GT. Numerical biharmonic analysis and some applications. International Journal of Solids and Structures 1967;3:309–32. [2] Jaswon MA, Maiti M. An integral formulation of plate bending problems. Journal of Engineering Mathematics 1968;2:83–93. [3] Bezine GP. Boundary integral formulation for plate flexure with arbitrary boundary conditions. Mechanics Research Communications 1978;5(4): 197–206. [4] Tottenhan H. The boundary element method for plates and shells. In: Banerjee PK, Butterfield R., editors. Developments in boundary element methods, 1979, p. 173–205. [5] Stern MA. A general boundary integral formulation for the numerical solution of plate bending problems. International Journal of Solids and Structures 1979;15:769–82. [6] Van der Weee¨n F. Application of the boundary integral equation method to Reissner’s plate model. International Journal for Numerical Methods in Engineering 1982;18:1–10. [7] Barcellos CA, Silva LHM. A boundary element formulation for Mindlin’s plate model. In: Brebbia CA, Venturini WS, editors. BETECH-87. Southampton: CMP; 1897. p. 123–30. [8] Aliabadi MH. Plate bending analysis with boundary element. Southampton: Computational Mechanics Publications; 1998. [9] Moshaiov A, Vorus WS. Elastoplastic plate bending analysis by a boundary element method with initial plastic moments. International Journal of Solids and Structures 1986;22(11):1213–29. [10] Chueiri LHM, Venturini WS. Elastoplastic BEM to model concrete slabs. In: Brebbia CA, et al., editors. Boundary elements XVII. CMP; 1995. p. 149–56.
[11] Providakis CP. Quasi-static analysis of viscoplastic plates by the boundaryelement method. Engineering Analysis with Boundary Elements 1993;11(4): 265–8. [12] Morjaria M. Inelastic analysis of transverse deflection of plates by the boundary element method. Journal of Applied Mechanics—Transactions of the ASME 1980;47:291. [13] Kamiya N, Sawaki Y. An integral equation approach to finite deflection of elastic plates. International Journal of Non-Linear Mechanics 1982;17(3): 187–94. [14] Ye T-Q, Liu Y-J. Finite deflection analysis of elastic plate by the boundary element method. Applied Mathematical Modelling 1985;9:183–8. [15] Tanaka M, Matsumoto T, Zheng Z-D. Incremental analysis of finite deflection of elastic plates via boundary-domain-element method. Engineering Analysis with Boundary Elements 1996;17:123–31. [16] Wang W, Ji X, Tanaka M. A dual boundary element approach for the problems of large deflection of thin elastic plates. Computational Mechanics 2000; 26:58–65. [17] Wen PH, Aliabadi MH, Young A. Large deflection analysis of Reissner plate by boundary element method. Computers & Structures 2005;83(10–11):870–9. [18] Purbolaksono J, Aliabadi MH. Large deformation of shear-deformable plates by the boundary element method. Journal of Engineering Mathematics 2005;51(3):211–30. [19] Dirgantara T, Aliabadi MH. A boundary element formulation for geometrically nonlinear analysis of shear deformable shells. Computer Methods in Applied Mechanics and Engineering 2006;195(33–36):4635–54. [20] Supriyono, Aliabadi MH. Boundary element method for shear deformable plates with combined geometric and material nonlinearities. Engineering Analysis with Boundary Elements 2006;30(1):31–42. [21] Fung YC. Foundations of solid mechanics. Englewood Cliff, NJ: Prentice-Hall, Inc.; 1964. [22] Brebbia CA, Telles JCF, Wrobel LC. Boundary element techniques. Theory and applications in engineering. Berlin and New York: Springer-Verlag; 1984. [23] Bui HD. Some remarks about the formation of three-dimensional thermoelastoplastic problems by integral equations. International Journal of Solids and Structures 1978;14:935–9. [24] Venturini WS, Paiva JB. Boundary elements for plate bending analysis. Engineering Analysis with Boundary Elements 1993;11:1–8. [25] Paiva JB, Venturini WS. Alternative technique for the solution of plate bending problems using the boundary element method. Advances in Engineering Software 1992;14(4):265–71. [26] Fernandes GR, Chaves EWV, Venturini WS. Plate bending boundary element formulation considering variable thickness. Engineering Analysis with Boundary Elements 1999;23:405–18. [27] Leite LGS, Venturini WS. Stiff and soft thin inclusions in two-dimensional solids by the boundary element method. Engineering Analysis with Boundary Elements 2003;29(3):257–67. [28] Leite LGS, Venturini WS. Stiff and soft thin inclusions in two-dimensional solids by the boundary element method. Engineering Analysis with Boundary Elements 2005;29(3):257–67. [29] Botta AS, Venturini WS, Benallal A. BEM applied to damage models emphasizing localization and associated regularization techniques. Engineering Analysis with Boundary Elements 2005;29:814–27. [30] Fernandes GR, Venturini WS. Building floor analysis by the boundary element method. Computational Mechanics 2005;35:277–91. [31] Fernandes GR, Venturini WS. Non-linear boundary element analysis of floor slabs reinforced with rectangular beams. Engineering Analysis with Boundary Elements 2007;31:721–37. [32] Bonnet M, Mukherjee S. Implicit BEM formulation for usual and sensitivity problems in elasto-plasticity using the consistent tangent operator concept. International Journal of Solids and Structures 1996;33:4461–80. [33] Benallal A, Fudoli CA, Venturini WS. An implicit BEM formulation for Gradient plasticity and localization phenomena. International Journal for Numerical Methods in Engineering 2002;53:1853–69. [34] Simo JC, Hughes TJR. Computational inelasticity. New York: Springer; 1998.