International Journal of Solids and Structures 44 (2007) 1038–1059

[Link]/locate/ijsolstr

Analysis of shear deformable plates with combined

geometric and material nonlinearities by boundary
element method
Supriyono, M.H. Aliabadi *

Department of Aeronautics, Faculty of Engineering, Imperial College London, Prince Consort Road,
London SW7 2BY, United Kingdom

Received 3 January 2006; received in revised form 3 April 2006

Available online 13 June 2006

Abstract

In this paper a boundary element formulation for analysis of shear deformable plates with combined geometric and
material nonlinearities by boundary element method is presented. The dual reciprocity method is used in dealing with
the geometric nonlinearity and domain discretization is implemented in dealing with material nonlinearity. The material
is assumed to undergo large deflection with small strains. The von Mises criteria is used to evaluate the plastic zone
and an elastic perfectly plastic material behaviour is assumed. An initial stress formulation is used to formulate the bound-
ary integral equations. A total incremental method is applied to solve the nonlinear boundary integral equations. Numer-
ical examples are presented to demonstrate the validity and the accuracy of the proposed method.
 2006 Elsevier Ltd. All rights reserved.

Keywords: Large deflection; Geometrical nonlinearity; Material nonlinearity; Plasticity; Dual reciprocity method

1. Introduction

The numerical analysis of plate type structures by mean of finite element method accounting from geomet-
rical and material nonlinearity has been developed significantly during the last two decades. The boundary
element method after the FEM is the second most utilized computational tool available in solid mechanics.
With increasing reliance on numerical method in place of experimental measurements, it is therefore essential
to develop alternative solution methods to FEM for the sake of validation.
The range of applications for boundary element method (BEM) in solid mechanics is getting broader (Alia-
badi, 2002). Recent developments have extended the application of the method to nonlinear problems in plate
bending analysis such as geometric nonlinearity, material nonlinearity and combined geometric and material
nonlinearities. In the formulation of the integral equation the nonlinear problems can normally be included in

*
Corresponding author.
E-mail address: [Link]@[Link] (M.H. Aliabadi).

0020-7683/$ - see front matter  2006 Elsevier Ltd. All rights reserved.
doi:10.1016/[Link].2006.06.004
 Supriyono, M.H. Aliabadi / International Journal of Solids and Structures 44 (2007) 1038–1059 1039

a domain integral, therefore two kinds of integrals appear in nonlinear problems, i.e., boundary integral and
domain integral. There are two main techniques for the evaluation of the domain integrals appearing in the
integral equation. The first one is a domain discretization method in which the domain is discretized into inter-
nal cells, so that the advantage of BEM, that is, a possibility of reduced dimensionality of the problem, dis-
appears. The second one is transformation the domain integral into equivalent boundary integrals so the
domain discretization is not required and the computational efficiency of BEM is maintained.
Domain discretization for dealing with domain integral of plate bending analysis in BEM has been reported
by many researchers. Lei et al. (1990) formulated an integral equation for geometrically nonlinear behaviour
for Reissner plate. To evaluate the domain integral appearing in the formulation, Lei et al. (1990) discretized
the domain into constant triangular cells. A similar work has been reported by Purbolaksono and Aliabadi
(2005), however they implemented constant rectangular cells instead of constant triangular cells. Their work
(Purbolaksono and Aliabadi, 2005) was focused on studying four methods of solution for the nonlinear prob-
lem which included total increment method, cumulative load incremental method, Euler method and nonlin-
ear system method. Purbolaksono and Aliabadi (2005) found the most efficient approach is the total increment
method proposed by Wen et al. (2005). Baiz and Aliabadi (2005) made an improvement in evaluating the effect
of geometric nonlinearity on the plate bending analysis by implementing 9-node quadrilateral cell and intro-
ducing free term factor in the integral equation.
Karam and Telles (1998) and later Ribeiro and Venturini (1998) used the domain discretization in dealing
with nonlinear due to material nonlinearity. Supriyono and Aliabadi (2006) were the first to formulate a BEM
for combined geometric and material nonlinearities for shear deformable plates. In dealing with the domain
integrals appearing in the formulation, Supriyono and Aliabadi (2006) discretized the domain into cells using
9-node quadrilateral cell. To solve the nonlinear system of equation the total incremental method as proposed
by Wen et al. (2005) was used. The formulation developed (Supriyono and Aliabadi, 2006) allows for large
deflection and small strain and elastic perfectly plastic material behaviour.
The dual reciprocity method (DRM) which transform domain integral into equivalent boundary integral
was introduced by Brebbia and Nardini (1983), for analyzing 2D elastodynamic problem. In plate bending
analysis the application of the transformation method of domain integral into equivalent boundary integral
can be traced through the works by Kamiya and Sawaki (1988), Silva and Venturini, 1985, Sawaki et al.
(1990), Wen et al. (2005). Kamiya and Sawaki (1988) analysed the bending problem of plate on Wrinkler-type
elastic foundation as a novel extension of the dual reciprocity method (DRM) to the problems concerning the
fourth order differential equation. Silva and Venturini (1985) applied DRM to a similar problem taking into
consideration the nonlinear effect of elastic foundation. Sawaki et al. (1990) analysed nonlinear bending of
thin plate based on the Berger equation. The domain integral which appears in the formulation was trans-
formed to the boundary using DRM and treated through a successive iteration scheme.
Wen et al. (2005) formulated the large deflection effect on the Reissner plate based on the general nonlinear
differential equations for large deflection. The nonlinear terms were treated as body forces and determinated by
approximation function. The domain integral was transformed to the boundary using DRM. Another contri-
bution of application of DRM in large deflection plate analysis can be found in the work by Wang et al. (1950).
In this paper the boundary integral equation for shear deformable plate with combined geometric and mate-
rial nonlinearities is presented. The two techniques for evaluating domain integral are employed. The dual rec-
iprocity method is used to evaluate the domain integral due to large deflection effect and domain discretization
technique is implemented in dealing with the nonlinear term due to plasticity. The nonlinear terms due to large
deflection were determinated by utilizing radial basis function and approximating the derivatives. The formu-
lation is based on an assumption that the material experiences large deflection and small strain (Stricklin et al.,
1972) and in the numerical implementation elastic perfectly plastic material is used. In order to solve the non-
linear system of equation the total incremental method is implemented as proposed by Wen et al. (2005).

2. Governing equations

In order to define a general formulation for combined geometrical and material nonlinearities of plate
bending, it is considered that plastic strains are only due to bending and membrane, hence total strain rates
can be defined as
1040 Supriyono, M.H. Aliabadi / International Journal of Solids and Structures 44 (2007) 1038–1059

v_ ab ¼ v_ eab þ v_ aab ð1Þ

e_ ab ¼ e_ eab þ e_ aab ð2Þ

and
c_ a3 ¼ c_ ea3 ð3Þ
where v_ ab are the total bending strain rates, e_ ab are the total in-plane strain rates, and c_ a3 are the shear strain
rates, respectively. The total bending strain rates consist of linear parts ðv_ eab Þ and inelastic parts ðv_ aab Þ. Similarly
total in-plane strain rates consist of linear parts ð_eeab Þ and nonlinear parts ð_eaab Þ. The nonlinear parts of (1) and
(2) are due to large deflection (nl) and plasticity (p), and they can be expressed as

v_ aab ¼ v_ nl _ pab
ab þ v ð4Þ
and
e_ aab ¼ e_ nl _ pab
ab þ e ð5Þ
The generalized strain rates above can be expressed in terms of generalized displacement rates as
w_ a;b þ w_ b;a
v_ ab ¼ ð6Þ
2
ðu_ a;b þ w_ a;b Þ þ w_ 3;b w3;a
e_ ab ¼ ð7Þ
2
and
c_ a3 ¼ w_ a þ w_ 3;a ð8Þ
where w_ a are rotation, u_ a are translation and w3 is deflection.
Naghdi (1956) derived the relationships between stress resultants and strains by using the Reissner’s vari-
ational theorem of elasticity (Reissner, 1950). The relationship can be written as
 
M_ ab ¼ D 1  m 2v_ ab þ 2m v_ cc dab  M _p ð9Þ
ab
2 1m
 
1m 2m
N_ ab ¼ B e_ ab þ e_ ba þ e_ cc dab  N_ pab ð10Þ
2 1m
and
Q_ a ¼ C c_ a3 ð11Þ
Eh3 Eh Ekh 5
where D ¼ ;B ¼1m2
and C ¼ 1m2
, 2ð1þmÞ
;k ¼ 6
.
The equilibrium equations are
M_ ab  Q_ a ¼ 0 ð12Þ
Q_ a;a þ ðN_ ab w_ 3;b Þ;a þ q_ 3 ¼ 0 ð13Þ
and
N_ ab;b ¼ 0 ð14Þ
where M _ ab ; Q_ a ; N_ ab are the moment stress resultants, the shear stress resultants, and membrane stress resul-
tants, respectively. q_ 3 is uniform load per unit area in the third direction.

3. Displacement integral equations

The governing equations together with the Betti’s reciprocal theorem are used to derive the displacement
integral equations. By considering equilibrium equations (12)–(14), it is possible to write the following
relationships:
 Supriyono, M.H. Aliabadi / International Journal of Solids and Structures 44 (2007) 1038–1059 1041
Z h i
_ ab;b  Q_ a ÞW  þ ðQ_ a;a þ ðN_ ab w_ 3;b Þ þ q_ 3 Þ þ W  dX ¼ 0
ðM ð15Þ
a ;a 3
X

and
Z
N_ ab;b U a dX ¼ 0 ð16Þ
X

where W i and U a are weighting functions and X is the domain where the integration is carried out. Integrat-
ing (15) by part and considering that P_ a ¼ M _ ab nb and p_ 3 ¼ Q_ a na and also considering relationships in (9)
gives
Z Z Z Z Z
p_ j W j dC þ ðN_ ab w_ 3;b Þ;a W 3 dX þ q_ 3 W 3 dX ¼ w_ j P j dC  ½ðM _   Q_  Þw_ a þ Q_  w_ 3  dX
ab;b a a;a
C X X C X
Z
 M _ p W  dX ð17Þ
ab a;b
X

By taking the state (Æ)* to correspond to concentrated generalized loads represented by Dirac delta function
D(X  X 0 ), one can obtain displacement integral equations for rotations (wa) and deflection (w3) from Eq.
(17) as follows:
Z Z Z
0 0

w_ i ðX Þ þ P ij ðX ; xÞw_ j ðxÞ dC ¼ W ij ðX ; xÞp_ j ðxÞ dC þ W i3 ðX0 ; XÞðN_ ab w_ 3;b Þ;a ðXÞ dX
 0

C C X
Z Z
þ W i3 ðX ; XÞq_ 3 ðXÞ dX þ viab ðX0 ; XÞM
 0 _ p ðXÞ dX ð18Þ
ab
X X

where x 2 C and X 2 X are field points on the boundary and the domain, respectively, W ij ðX0 ; xÞ; P ij ðX0 ; xÞ,
and viab ðX0 ; XÞ are the fundamental solutions of displacements, tractions, and strains, respectively. The expres-
sions of the kernels are given in Appendix A.
The in-plane displacement integral equations can be derived in a similar way. Integrating Eqs. (16) by part
gives
Z Z
N_ ab nb U a dC  N_ ab U a;b dX ¼ 0 ð19Þ
C X

Utilizing the relationships in (10), the integral equation (19) can be written as
Z Z Z Z
_ta U a dC ¼ N ab U a;b dX þ N ab U a;b dX  N_ pab U a;b dX
_ e  _ nl 
ð20Þ
C X X X

By taking the state (Æ)* to correspond to concentrated generalized loads represented by Dirac delta function
D(X  X 0 ), the displacement integral equations for in-plane displacement can be stated as
Z Z Z
u_ h ðX0 Þ þ T ha ðX0 ; xÞu_ a ðxÞ dC ¼ U ha ðX0 ; xÞt_a ðxÞ dC þ U ha ðX0 ; XÞN_ nl
ab;b ðXÞ dX
C C X
Z
þ ehab ðX0 ; XÞN_ pab ðXÞ dX ð21Þ
X

where U ha ðX0 ; xÞ; T ha ðX0 ; xÞ,

and ehab ðX0 ; XÞ
are the fundamental solutions of displacements, tractions, and
strains, respectively. The expressions of the kernels are given in Appendix A.
By taking the point X 0 to the boundary, that is X 0 ! x 0 2 C, assuming that displacements uj satisfy Holder
continuity, Eq. (18) can be written as follows:
Z Z Z
0  0
cij wi ðx Þ þ - P ij ðx ; xÞw_ j ðxÞ dC ¼ W ij ðx ; xÞp_ j ðxÞ dC þ W i3 ðx0 ; XÞðN_ ab w_ 3;b Þ;a ðXÞ dX
 0

C C X
Z Z
þ W i3 ðx ; XÞq_ 3 ðXÞ dX þ viab ðx0 ; XÞM
 0 _ p ðXÞ dX ð22Þ
ab
X X
1042 Supriyono, M.H. Aliabadi / International Journal of Solids and Structures 44 (2007) 1038–1059

and Eq. (21) can be written as

Z Z Z
cha u_ h ðx0 Þ þ - T ha ðx0 ; xÞu_ a ðxÞ dC ¼ U ha ðx0 ; xÞt_a ðxÞ dC þ U ha ðx0 ; XÞN_ nl
ab;b ðXÞ dX
C C X
Z
þ ehab ðx0 ; XÞN_ pab ðXÞ dX ð23Þ
X
R
where - denotes a Cauchy principal value integral. Eqs. (22) and (23) constitute the boundary displacement
integral equations for nonlinear plate bending problem.

4. Stress integral equations

The stress resultants at the domain point X 0 can be evaluated by differentiating Eqs. (18) and (21) with
respect to the source point X 0 and then substituting the resulting expression into Eqs. (9)–(11). The differen-
tiation can be applied directly to the fundamental solutions of Eqs. (18) and (21), except for viab ðX0 ; XÞ and
ehab ðX0 ; XÞ, for which the Leibnitz formula must be considered. The stress integral equations for moment stress
resultants can be stated as
Z Z Z
M_ ab ðX0 Þ ¼ W  ðX0 ; xÞp_ k ðxÞdC  P  ðX0 ; xÞw_ k ðxÞdC þ W  ðX0 ; X Þq_ 3 ðXÞdX
abk abk ab3
C C X
Z Z _p _ p dab 
0 0 p ½2ð1 þ mÞM ab þ ð1  3vÞM hh
þ  _ _ 3;b Þ ðXÞdX þ  _
W ab3 ðX ; XÞðN cb w ;c vabch ðX ; XÞM ch ðXÞdX 
X X 8
ð24Þ
and shear stress resultants can be written as
Z Z Z
Q_ a ðX0 Þ ¼ W 3bk ðX0 ; xÞp_ k ðxÞ dC  P 3bk ðX0 ; xÞw_ k ðxÞ dC þ W 3b3 ðX0 ; XÞq_ 3 ðXÞ dX
C C X
Z Z
þ W 3b3 ðX0 ; XÞðN_ cb w_ 3;b Þ;c ðXÞ dX þ v3bch ðX0 ; XÞM _ p ðXÞ dX ð25Þ
ch
X X
finally membrane stress resultants can be expressed as
Z Z Z
N_ ab ðX0 Þ ¼ U abc ðX0 ; xÞt_c ðxÞ dC  T abc ðX0 ; xÞu_ c ðxÞ dC þ U abc ðX0 ; XÞN_ nl
ch;h ðXÞ dX
C C X
Z
½2ð1 þ mÞN_ pab þ ð1  3vÞN_ phh dab 
þ eabch ðX0 ; XÞN_ pch ðXÞ dX  ð26Þ
X 8
where the kernels W ibk ; P ibk are linear combination of the first derivatives of W ij and P ij . The kernels U abc ; T abc
are the linear combination of the first derivatives of U ab and T ab . The kernels vibch ; eibch are the linear combi-
nation of the first derivatives of viab and eabc . The free terms appear in Eqs. (24) and (26) arising from using
Leibnitz formula. The expressions of the kernels are listed in Appendix B.

5. Transformation of domain integrals and determination of nonlinear terms due to large deflection

5.1. Domain integrals

In displacement and stress integral equation above, there are four domain integrals due to large deflection
effect as follows:
Z
I1 ¼ W i3 ðx0 ; XÞðN_ ab w_ 3;b Þ;a ðXÞ dX
X
Z
I2 ¼ W ib3 ðX0 ; XÞðN_ cb w_ 3;b Þ;c ðXÞ dX
ZX
I3 ¼ U ha ðx0 ; XÞN_ nl
ab;b ðXÞ dX
X
 Supriyono, M.H. Aliabadi / International Journal of Solids and Structures 44 (2007) 1038–1059 1043

and
Z
I4 ¼ U abc ðX0 ; XÞN_ nl
ch;h ðXÞ dX ð27Þ
X

The dual reciprocity method (DRM) is used to transform the domain integrals into equivalent boundary
integrals. The transformation is made after approximating the nonlinear terms by a set of radial basis function
Fm(rm) associated with M points Xm distributed within domain X. Thus, the nonlinear term ðN_ ab w_ 3;b Þ;a in
domain integrals I1 and I2 are approximated by
X
M
ðN_ ab w_ 3;b Þ;a ¼ F m ðrm Þam ð28Þ
m¼1

where a are a set of unknown coefficients and r = jX  Xmj is the distance between the field point X and the
m

interior point Xm. If the nonlinear term ðN_ ab w_ 3;b Þ;a is written as a vector f, Eq. (28) can be inverted to obtain
am, i.e.,
a ¼ F1 f ð29Þ
If the particular solution w
^ ml satisfy the differential equation
^ mkl ¼ f ðrÞdil
Lik w ð30Þ
the domain integrals I1 and I2 can be transformed into equivalent boundary integrals as
X
M  Z Z 
m 0 m 0 m 0 m
I1 ¼ ak cij ðx Þ^
wjk  W ij ðx ; xÞ^
pjk ðxÞ dC þ - P ij ðx ; xÞ^
wjk ðxÞ dC ð31Þ
m¼1 C C

and
X
M  Z Z 
I2 ¼ amk b ijk ðx0 Þ 
M 0
W ijl ðx pmlk ðxÞ dC
; xÞ^ 0
þ - P ijl ðx ; xÞ^ m
wlk ðxÞ dC ð32Þ
m¼1 C C

b 3a ¼ Q
where shear force particular solution M b a . The particular solution w ^ mjk and the associated traction ^pmlk and
b
M ijk are given in Appendix C.
In the similar way, the domain integrals I3 and I4 can be transformed into equivalent boundary integral as
XM  Z Z 
m 0 m 0 m 0 m
I3 ¼ wbg  U ab ðx ; xÞ^tbg ðxÞ dC þ - T ab ðx ; xÞ^ubg ðxÞ dC
ag cab ðx Þ^ ð33Þ
m¼1 C C

and
X
M  Z Z 
I4 ¼ amg N^ ab ðx0 Þ  U abc ðx0 ; xÞ^tm ðxÞ dC þ - T abc ðx0 ; xÞ^um ðxÞ dC ð34Þ
cg cg
m¼1 C C

^mab and the associated traction tmab are given in Appendix C.

The particular solution u

5.2. Nonlinear terms

The nonlinear term in the above equations are N_ nl _ _ 3;b Þ. The term N_ ab can be determined from Eq.
ab and ðN ab w
_ nl
(26), whereas N ab and w_ 3;b are calculated as follows:
1m  m 
N_ nl
ab ¼ B w_ 3;b w_ 3;a þ w_ 3;c w_ 3;c dab ð35Þ
2 1m
and
Z Z Z
w_ 3;c ðX0 Þ þ P 3j;c ðX0 ; xÞw_ j ðxÞ dC ¼ W 3j;c ðX0 ; xÞp_ j ðxÞ dC þ W 33;c ðX0 ; XÞðN_ ab w_ 3;b Þ;a ðXÞ dX
C C X
Z
 0
þ W 33;c ðX ; XÞq_ 3 ðXÞ dX ð36Þ
X
1044 Supriyono, M.H. Aliabadi / International Journal of Solids and Structures 44 (2007) 1038–1059

To calculate the derivatives of the nonlinear terms, the nonlinear terms are approximated by
X
M
Sðx1 ; x2 Þ ¼ f ðrÞam ð37Þ
m¼1

where S can be either N_ nl _ _ 3;b Þ; am are coefficient which can be obtained by the values at selected point
ab or ðN ab w
and approximating function
pffiffiffiffiffiffiffiffiffiffiffiffiffiffi
f ðrÞ ¼ c2 þ r2 ð38Þ
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
2 2
where c is a constant and r ¼ ðx1  xm1 Þ þ ðx1  xm1 Þ . Then, the derivatives of the nonlinear terms are cal-
culated from the formula
XM
ðxa  xma Þam
S ;a ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð39Þ
m¼1 c2 þ r 2

6. Elastoplastic constitutive equations

In elastoplastic analysis, the existence of a yield function U is admitted. The yield function U is expressed in
terms of the stresses rab and in term of a hardening parameter k. During the loading that produces yielding,
the stresses rab must remain at the yield surface, so that the following equation is satisfied:
Uðrab ; kÞ ¼ f ðrab Þ  WðkÞ ¼ re  r0 ð40Þ
where re is the equivalent stress calculated using von Mises yield criteria in this work and r0 is the uniaxial
yield stress. In the case where the membrane and moment stresses exist at the same time on the analysis of
Reissner plate, the equivalent stress re and the uniaxial stress r0 can be stated as
1 4
re ¼ N e þ 2 M e ð41Þ
h h
and
1 4
r0 ¼ N 0 þ 2 M 0 ð42Þ
h h
where Ne and Me are the equivalent membrane and moment stress, respectively. N0 and M0 are the uniaxial
membrane and moment stress, respectively. It is considered here that whenever the equivalent stress re at any
point reaches the yield stress r0, the whole cross section yields simultaneously.
For von Mises yield criteria, Ne and Me are calculated as
pffiffiffiffiffiffiffi
N e ¼ 3J 2 ð43Þ
and
qffiffiffiffiffiffiffi
Me ¼ 3J 02 ð44Þ

J2 and J 02 are analogous to the second invariant of the deviatoric stress that can be expressed as
J 2 ¼ 12N 0ij N 0ij ð45Þ
and
J 02 ¼ 12M 0ij M 0ij ð46Þ
with
N 0ij ¼ N ij  13dij N kk ð47Þ
and
M 0ij ¼ M ij  13dij M kk ð48Þ
 Supriyono, M.H. Aliabadi / International Journal of Solids and Structures 44 (2007) 1038–1059 1045

The stress–strain relationships for the behaviour after yielding can be stated in the incremental form as
1
dM ab ¼ C ep
abch dvch  C ablq alq agf dgf ð49Þ
c0
and
1
dN ab ¼ C ep
abch dech  C ablq alq agf dgf ð50Þ
c0
where
1
C ep
abch ¼ C abch  C ablq alq agf C gfch ð51Þ
c0
where Cabch represents the components of fourth order isotropic tensor of elastic constants, c 0 and ach can be
written as
c0 ¼ aab C abch ach þ H 0 ð52Þ
where
oU
ach ¼
oM ch
for moment and for membrane ach can be stated as
oU
ach ¼ ð53Þ
oN ch
where H 0 ¼ ov
oW 0
p for moment and H ¼
oW
oe
0
p membrane. H is called the slope of the stress–plastic strain curve.
e e

7. Discretization and system of equations

Substituting Eqs. (31) and (32) into Eqs. (22) and (23), one gets
Z
cij_wi ðx Þ þ - P ij ðx0 ; xÞw_ j ðxÞ dC
0

C
Z X
M  Z Z 
¼ W ij ðx0 ; xÞp_ j ðxÞ dC þ amk cij ðx0 Þ^
wmjk  W ij ðx0 ; xÞ^pmjk ðxÞ dC þ - P ij ðx0 ; xÞ^
wmjk ðxÞ dC
C m¼1 C C
Z Z
þ W i3 ðx0 ; XÞq_ 3 ðXÞ dX þ _ p ðXÞ dX
viab ðx0 ; XÞM ð54Þ
ab
X X

and
Z
c_ha uh ðx0 Þ þ - T ha ðx0 ; xÞu_ a ðxÞ dC
C
Z X
M  Z Z 
¼ U ha ðx0 ; xÞt_a ðxÞ dC þ amg cab ðx0 Þ^
wmbg  U ab ðx0 ; xÞ^tmbg ðxÞ dC þ - T ab ðx0 ; xÞ^umbg ðxÞ dC
C m¼1 C C
Z
þ ehab ðx0 ; XÞN_ pab ðXÞ dX ð55Þ
C

In order to analyse the problem by the boundary element method, the integral Eqs. (54) and (55) are dis-
cretized. The boundary C is divided into boundary elements and the domain X where the existence of plastic
deformation is expected is divided into cells. The boundary is discretized using quadratic isoparametric ele-
ments both continuous and semi-discontinuous elements are used and the domain is discretized using 9-node
quadrilateral (continuous and semi-discontinuous) cells. The semi-discontinuous boundary elements are used
to avoid corner problems in the boundary and the semi-discontinuous cells to avoid the coincident nodes
1046 Supriyono, M.H. Aliabadi / International Journal of Solids and Structures 44 (2007) 1038–1059

between boundary and cell. The internal values are taken at the cell nodes which will avoid calculation of
deflection derivative on the boundary nodes. The deflection derivative is used to calculate the nonlinear terms
due to large deflection. In dealing with the large deflection effect, the DRM points are taken the same as the
cell points.
After discretization and point collocation on the boundary as well as in the domain, Eqs. (54) and (55) can
be written in the matrix form as
 w 
 w 
( _ ) " _ w #  w ( p )
H 0 w_ G 0 p_ b f T 0 M_
u ¼ u þ þ u þ u ð56Þ
0 H u_ 0 G t_ 0 f_ 0 T N_ p

where [H] and [G] are the well-known boundary element influence matrices, [T] is the influence matrices for
plasticity. The superscripts w and u show the plate and the in-plane modes, respectively, fwg; _ fug; _ ft_g
_ fpg;
_ _
are the displacement and the traction rate vectors. fbg and ffg are the load rate vectors due to external load
and large deflection, respectively. fM _ p g; fN_ p g are the bending and membrane stress resultant term, respec-
tively. After imposing boundary condition, Eq. (56) can be written as
 w ( p )
T 0 M_
½Afxg _
_ ¼ ffg þ ð57Þ
u
0 T N_ p

where [A] is the system matrix, fxg

_ is the unknown vector and ffg_ is the vector of prescribed boundary values.
Analogously, the stress integral equations can be presented in matrix form as
8 9 2 wa 3 2 wa 3 8 9 2 wa 3 2 wa 3
> _ 0
0
> _a f_ T þ Ewa 0 ( )
<M > = H
w_
G
p_ <b > = M_p
6 w3 7 6 w3 7 6 7 6 w3 7
Q_ ¼ 4 H 0 5 þ4G 0 5 þ b_ 3 þ 4 f 5 þ 4 T
_
w3
0 5
: _ >
> ; u u_ u t_ : >
> ; N_ p
N 0 H 0 G 0 f u_
0 T þ Eu
u

ð58Þ
_ fQg;
where fMg; _ fN_ g are the vectors of bending stress resultants, shear stress resultants and membrane stress
resultants, respectively. Superscripts wa and w3 denote the bending and shear modes, respectively.

8. Solution algorithm

To achieve a linearization of the nonlinear integral equations, the total increment method as proposed by
Wen et al. (2005) is applied. For the first step of the incremental procedure, the nonlinear terms due to both
large deflection and plasticity are set equal to zero, then for the step of kth the approximations of the nonlinear
terms are estimated using the (k  1)th step, so
h i h i
ðN_ cb w_ 3;b Þ;a ¼ ðN_ cb w_ 3;b Þ;a
k k1
h i h i
ðN_ cb Þ;a ¼ ðN_ cb Þ;a
nl nl
k k1
_ p Þ ¼ ðM
ðM _ pÞ
k k1

and
ðN_ p Þk ¼ ðN_ p Þk1
To evaluate the plastic zone of the model, the von Misses criteria is used.
By considering that
_e ¼M
M _p
_ ab þ M
ab ab

and
N_ eab ¼ N_ ab þ N_ pab
 Supriyono, M.H. Aliabadi / International Journal of Solids and Structures 44 (2007) 1038–1059 1047

Eq. (58) can be written as

8 _ e9 2 wa 3 2 wa 3 8 9 2 3 2 wa 3
H 0 ( ) G 0 ( ) > b_ a > _
f wa T þ Ewa þ I 0 ( )
<M >
> = >
< >
= 6 w3 7 _p
6 7 w_ 6 w3 7 p_ 6 7 M
Q_ e ¼ 4 H
w3
0 5 þ 4G 0 5 þ b_ 3 þ 6 _ 7
4f 5 þ 4 Tw3 0 5
>
: e> ; u_ _
t >
> >
> N_ p
u u : ; u u
N_ 0 H 0 G 0 f u_
0 T þE þI
ð59Þ

where M _ e and N_ e are elastic moment stress resultants and elastic membrane stress resultant, respectively, I is
an identity matrix.
After computation of all the matrices and known vectors, every load step the following system matrices are
solved:
 w ( p )
T 0 M_ þ DM _p
½Afxg_ ¼ ffg_ þ ð60Þ
0 Tu N_ p þ DN_ p
and
8 9 2 wa 3 2 wa 3 8 9 2 wa 3 2 wa 3
> _e 0
0
> _a f_ T þ Ewa þ I 0 ( )
<M > = H
w_
G
p_ <b > = M_ p þ DM_p
6 w3 7 6 w3 7 6 7 6 w3 7
Q_ e ¼ 4 H 0 5 þ 4G 0 5 þ b_ 3 w3
_
þ 4f 5 þ 4 T 0 5
: _e>
> ; u u_ u t_ : >
> ; u_ u u N_ p þ DN_ p
N 0 H 0 G 0 f 0 T þ E þ I
ð61Þ

where DM_ p and DN_ p denote the increment plastic resultants.

The nonlinear terms due to plasticity, can be calculated as shown by Karam and Telles (1998). Assuming an
incremental fictitious ‘‘elastic moment and membrane’’ we have

Fig. 1. Square plate mesh: (a) BEM mesh and (b) FEM mesh.

Fig. 2. Circular plate mesh: (a) BEM mesh and (b) FEM mesh.
1048 Supriyono, M.H. Aliabadi / International Journal of Solids and Structures 44 (2007) 1038–1059

dM eab ¼ C abch dvch ð62Þ

and
dN eab ¼ C abch dech ð63Þ
Taking into account Eqs. (62) and (63), Eqs. (49) and (51) can be written as
1
dM ab ¼ dM eab  C ablq alq agf dM egf ð64Þ
c0
and
1
dN ab ¼ dN eab  C ablq alq agf dN egf ð65Þ
c0

7
Increment (ΔQ)= 0.144
6

4
Q

81 DRM 4x4 Cells

3
121 DRM 5x5 Cells

2 149 DRM 6x6 Cells

FEM 8x8 Cells

1

0
0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6

Fig. 3. Comparison of center point deflection of different BEM meshes to the FEM for DQ = 0.288.

8
149 DRM Points, 6x6 Cells
7

5
Q

4 Incr.=0.36
Incr.=0.24
3
Incr.=0.18
Incr.=0.144
2

0
0 0.5 1 1.5 2 2.5 3
W

Fig. 4. Comparison of center point deflection of different increment sizes for mesh C.
 Supriyono, M.H. Aliabadi / International Journal of Solids and Structures 44 (2007) 1038–1059 1049

The steps in calculating the nonlinear terms due to plasticity can be described as follows:

(i) Solve Eq. (61) to obtain the moment and membrane elastic stress increments.
(ii) Calculate the true moment and membrane stress increments by solving Eqs. (64) and (65).
(iii) Calculate the nonlinear terms due to plasticity, by
DM pab ¼ DM eab  DM ab ð66Þ

and
DN pab ¼ DN eab  DN ab ð67Þ
(iv) Accumulate the values of nonlinear terms due to plasticity:

5
Q

3 Comb.
Plastic
2 Linear
Large Defl.
1

0
0.0 0.5 1.0 1.5 2.0
W

Fig. 5. Comparison of center point deflection of different BEM analysis for square plate.

300 MPa
300MPa

300MPa

300 MPa

Fig. 6. The contour of plastic zone after plasticity takes place at Q = 3.168.
1050 Supriyono, M.H. Aliabadi / International Journal of Solids and Structures 44 (2007) 1038–1059

ðM pab Þk ¼ ðM ab Þk1 þ ðDM pab Þk ð68Þ

and
ðN pab Þk ¼ ðN ab Þk1 þ ðDN pab Þk

9. Numerical examples

Two numerical examples are presented to demonstrate the validity and the accuracy of the proposed BEM
formulation. The configurations studied are square and circular plates subjected to a uniform distributed load
q with all sides simply supported. The plate discretization are shown in Figs. 1 and 2. Comparisons are made
with numerical solutions obtained using the finite element method. In the following examples the parameters
qa4 wmax
Q ¼ Eh 4 and W ¼ h
, where E is the modulus elasticity of the material, wmax is the deflection of the center
point of the plates, h is the thickness of the plates and a is the side for square plate or radius for circular plate.
Example 1. Simply supported square plate subjected to a uniformly lateral distributed load q.
The square plate has properties of E = 200 GPa, m = 0.3 and ryield = rult = 300 MPa, where ryield, rult are
yield and ultimate stresses, respectively. Material is assumed to be an elastic perfectly plastic and ha ¼ 0:05.
Three different meshes are used to study the influence of the discretization on the results. The first mesh A
has 32 boundary elements of quadratic isoparametric element, 4 · 4 9-node cells and 81 DRM nodes. The
second mesh B has 32 quadratic isoparametric boundary elements, 5 · 5 9-node cells and 121 DRM nodes.
The third mesh C has 32 quadratic isoparametric boundary elements, 6 · 6 9-node cells and 149 DRM
nodes. For a selected mesh four different increment sizes (DQ) which are 0.36, 0.24, 0.18 and 0,144 are
implemented.
The deflections at the center point of the plate for every load step/increment are presented in Figs. 3–5.
Firstly, in Fig. 3 the results of combined nonlinear BEM analysis of the three meshes for DQ = 0.144 are com-
pared to the FEM results. The FEM results are obtained from ABAQUS on which the model is discretized into
8 · 8 elements of 8-node shell. It can be seen, BEM model is able to simulate the problem of combined geomet-
rical and material nonlinearities. Both BEM and FEM give similar load at the first yield when structure starts
to collapse, however after a certain load, the deflections obtained by BEM are larger than the ones obtained by
FEM. The increased number of domain cells and DRM points improves the accuracy of the results. Further-
more, it is shown in Fig. 4 that for the total incremental method which is implemented to deal with the non-
linear system of equation, the larger number of load increments (smaller DQ) also gives better results.

6.0

5.0

4.0

3.0 Incr.=0.112, 197 DRM, 44 cells

Q

Incr.=0.112, 145 DRM, 32 cells

Incr.=0.28, 197 DRM, 44 cells
2.0
FEM

1.0

0.0
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
W

Fig. 7. Comparison of center point deflection of different BEM meshes and increment size to the FEM results.
 Supriyono, M.H. Aliabadi / International Journal of Solids and Structures 44 (2007) 1038–1059 1051

The results of combined nonlinear BEM analysis of mesh C compared to the linear, large deflection and
plasticity results also obtained by BEM analysis are presented in Fig. 5. The results for combined case are
found to be closer to the elastoplastic results, hence the influence of the plasticity is greater than large deflec-
tion in this particular example. Finally, the plastic zone after plasticity takes place at Q = 3.168 is presented in
Fig. 6. It can be seen that plastic zone is started from the corner and the center of the plate.
Example 2. Simply supported circular plate subjected to a uniformly lateral distributed load q.
The circular plate has the same material properties as square plate in Example 1, the ratio between the
thickness and the radius ðhrÞ is 0.05. Two different meshes are implemented. The first one, A, has 16 boundary
elements, 32 cells and 145 DRM points and the second one, B, has 32 boundary elements, 44 cells and 197
DRM points. For a selected mesh two different increment sizes (DQ) which are 0.28 and 0.112 are
implemented.

6.0

5.0

4.0

3.0
Q

Comb.
2.0 Plastic
Large defl.
1.0 Linear

0.0
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
W

Fig. 8. Comparison of center point deflection of different BEM analysis for circular plate.

300 MPa

300MPa

Fig. 9. The contour of plastic zone after plasticity takes place at Q = 3.024.
1052 Supriyono, M.H. Aliabadi / International Journal of Solids and Structures 44 (2007) 1038–1059

The results are presented in Figs. 7–9. Fig. 7 presents comparison of center point deflection of various BEM
meshes and increment sizes to the FEM results. The comparison of various BEM analysis which are linear,
large deflection, combined nonlinearities and plasticity are presented in Fig. 8. Finally, Fig. 9 presents the
contour of plastic zone after plasticity takes place at Q = 3.024. From Figs. 7 and 8, it can be seen that the results
have similar trends as the square plate. The more number of domain cells, DRM point and increment the
better results can be obtained. For this case, the simply supported circular plate, the plastic zone is concen-
trated at the center of the plate.

10. Conclusions

The application of BEM to combined geometric and material nonlinearities for shear deformable plate bend-
ing analysis was presented where dual reciprocity method was used in dealing with large deflection effect and
domain cell technique was employed in dealing with plasticity. The formulations of the boundary integral equa-
tions were presented in detail including the integral equations used to calculate the internal values and transfor-
mation of the domain integrals due to large deflection to the equivalent boundary integrals. The numerical
implementation of the formulation was performed. From the results obtained, it can be concluded that

• the dual reciprocity method (DRM) and cell discretization technique can be coupled for dealing with com-
bined nonlinearities in order to increase the efficiency of BEM computation;
• good agreement with FEM results can be achieved with relatively coarser mesh;
• the total incremental method was shown to be an efficient approach for this problem as repeated solution of
system of equations is not required and the nonlinear terms are updated by back substitution.

Acknowledgement

The first author gratefully acknowledge the Muhammadiyah University of Surakarta and the management
and staffs of TPSDP project of Muhammadiyah University of Surakarta for their services and supports.

Appendix A

The expressions of the kernels in Eq. (18) are presented here:

1  
W ab ¼ ½8BðzÞ  ð1  mÞð2 ln z  1Þdab  ½8AðzÞ þ 2ð1  mÞr;a r;b
8pDð1  mÞ
1
W a3 ¼ W 3a ¼ ð2 ln z  1Þrr;a ðA:1Þ
8pD
1
W 33 ¼ ½ð1  mÞz2 ðln z  1Þ  8 ln z
8pDð1  mÞk2

and
1 
P ca ¼ ð4AðzÞ þ 2zK 1 ðzÞ þ 1  mÞ dac r;n þ r;a nc þ ð4AðzÞ þ 1 þ mÞr;c na
4pr 
2ð8AðzÞ þ 2zK 1 ðzÞ þ 1  mÞr;a r;c r;n
k2
P c3 ¼ ½BðzÞnc  AðZÞr;c r;n  ðA:2Þ
2p   
ð1  mÞ ð1 þ mÞ
P 3a ¼ 2 ln z  1 na þ 2r;a r;n
8p ð1  mÞ
1
P 33 ¼ r;n
2pr
 Supriyono, M.H. Aliabadi / International Journal of Solids and Structures 44 (2007) 1038–1059 1053

and
1
vcab ¼  ½ð8AðzÞ þ 4zK 1 ðzÞ þ 2  2mÞðdbc r;a þ dac r;b Þ
8pDð1  mÞr
 4ð8AðzÞ þ 2zK 1 ðzÞ þ 1  mÞr;a r;b r;c þ 2ð4AðzÞ þ 1  mÞdab r;c  ðA:3Þ
1
v3ab ¼  ½ð2 ln z  1Þdab þ 2r;a r;b 
8pD
where
 
2 1
AðzÞ ¼ K 0 ðzÞ þ K 1 ðzÞ 
z z
  ðA:4Þ
1 1
BðzÞ ¼ K 0 ðzÞ þ K 1 ðzÞ 
z z
in which K0(z) and K1(z) are the modified Bessel functions, z = kr, r is the absolute distance between the source
and the field points, r;a ¼ rra , where ra = xa(x)  xa(x 0 ) and r,n = r,ana.
Expanding the modified Bessel functions for small arguments:
h  
zi h z i ðz2 =4Þ
2
1 z ðz2 =4Þ
K 0 ðzÞ ¼ c  ln þ c þ 1  ln þ c þ 1 þ  ln
2 2 ð1!Þ2 2 2 ð2!Þ2
  2 3
1 1 z ðz =4Þ
þ c þ 1 þ þ  ln þ 
2 3 2 ð3!Þ2
    ðA:5Þ
1=2 3=2
1 1 z ðz2 =4Þ 1 z ðz2 =4Þ
K 1 ðzÞ ¼  c þ  ln  c þ 1 þ  ln
z 2 2 0!1! 4 2 ð1!2!Þ
  2 5=2
1 1 z ðz =4Þ
 c þ 1 þ þ  ln þ 
2 6 2 ð2!3!Þ
where c = 0.5772156649 is the Euler constant. Substituting Eqs. (A.5) into Eqs. (A.4) and taking the limit as
r ! 0:
1
lim AðzÞ ¼ 
r!0 2  ðA:6Þ
1 z 1
lim BðzÞ ¼  lim ln þ c þ
r!0 2 r!0 2 2
As can be seen that A(z) is a smooth function, whereas B(z) is a weakly singular O(ln r). Therefore, W ij are
weakly singular and P ij are strongly singular Oð1r Þ.
The expressions for the kernels in Eq. (21) are
 
1 1
U ha ¼ ð3  mÞ ln dha þ ð1 þ mÞr;a r;h ðA:7Þ
4pBð1  mÞ r
1
T ha ¼  fr;n ½ð1  mÞdha þ 2ð1 þ mÞr;a r;h  þ ð1  mÞ½nh r;a  na r;h g ðA:8Þ
4pr
1
ehab ¼ ½ð1  mÞðdha r;b þ dhb r;a Þ þ ð1 þ mÞðdab r;h þ 2r;a r;b r;h Þ ðA:9Þ
4pBð1  mÞr
U ha are weakly singular Oðln 1r Þ and T ha are strongly singular Oð1r Þ.

Appendix B

The calculation of the stresses at the internal points, one can use Eqs. (24)–(26). The kernel expressions of
the equations are given by
1054 Supriyono, M.H. Aliabadi / International Journal of Solids and Structures 44 (2007) 1038–1059

1  
W abc ¼ ð4AðzÞ þ 2zK 1 ðzÞ þ 1  mÞðdbc r;a þ dac r;b Þ2ð8AðzÞ þ 2K 1 ðzÞ þ 1  mÞr;a r;b r;c þ ð4AðzÞ þ 1 þ mÞdab r;c
4pr

and
  
ðm  1Þ ð1 þ mÞ
W ab3 ¼ 2 ln z  1 dab þ 2r;a r;b
8p 1m
k2 ðB:1Þ
w3bc ¼ ½BðzÞdcb  AðzÞr;b r;c 
2p
1
W 3b3 ¼ r;b
2pr
Dð1  mÞ 
P abc ¼ ð4AðzÞ þ 2K 1 ðzÞ þ 1  mÞðdca nb þ dcb na Þ þ ð4AðzÞ þ 1 þ 3mÞdab nc Þ  ð16AðzÞ þ 6zK 1 ðzÞ
4pr2
þ z2 K 0 ðzÞ þ 2  2mÞ  ½ðna r;b  nb r;a Þr;c þ ðdbc r;a þ dac r;b Þr;n   2ð8A9zÞ þ 2ðK 1 ðzÞ þ 1 þ mÞðdab r;c r;n
þ nc r;a r;b Þ þ 4ð24A9zÞ þ 8zK 1 ðzÞ þ z2 ðK 0 ðzÞ þ 2  2mÞr;a r;b r;c r;n

and

Dð1  mÞk2
P ab3 ¼ ½ð2AðzÞ þ zK 1 ðzÞÞðna r;b þ nb r;a Þ  2ð4AðzÞ þ zK 1 ðzÞÞr;a r;b r;n þ 2AðzÞdab r;n 
4pr
Dð1  mÞk2  
ðB:2Þ
P 3bc ¼ 2ðAðzÞ þ zK 1 ðzÞÞðdcb r;n þ r;c nb Þ þ 2AðzÞnc r;b  2ð4AðzÞ þ zK 1 ðzÞÞr;b r;c r;n
4pr
Dð1  mÞk2 2
P 3b3 ¼ ½ðz BðzÞ þ 1Þnb  ðz2 AðzÞ þ 2Þr;b r;n 
4pr2
r  
Qab ¼ ð4 ln z  3Þ½ð1  mÞðna r;b þ nb r;a Þ þ ð1 þ 3mÞdab r;n  þ 4½ð1  mÞr;a r;b þ mdab r;n
64p
and
1
Q3b ¼ ½ð2 ln z  1Þnb þ 2r;b r;n  ðB:3Þ
8p
1 
vabch ¼ 4ð4AðzÞ þ 2zK 1 ðzÞ þ 1  mÞðdca dhb þ dcb dha Þ þ 4ð4AðzÞ þ 1 þ mÞdab dch  4ð16AðzÞ þ 6zK 1 ðzÞ
16pr2
þ z2 K 0 ðzÞ þ 2  2mÞ½dha r;b r;c þ dhb r;a r;c þ dca r;b r;h þ dab r;a r;h   8ð8AðzÞ þ 2zK 1 ðzÞ þ 1 þ mÞdab r;c r;h

 8ð8AðzÞ þ 2zK 1 ðzÞ þ 1 þ mÞdch r;a r;b þ16ð24A9zÞ þ ð8zK 1 ðzÞ þ z2 K 0 ðzÞ þ 2  2mÞr;a r;b r;c r;h

and
k2  
v3bch ¼ ð8AðzÞ þ 4zK 1 ðzÞÞðdcb r;h þ dhb r;c Þ  4ð8AðzÞ þ 2zK 1 ðzÞÞr;b r;c r;h þ 8AðzÞdch r;b ðB:4Þ
16pr
1
U abc ¼  ½ð1  mÞðdca r;b þ dcb r;a  dab r;c Þ þ 2ð1 þ mÞr;a r;b r;c  ðB:5Þ
4pr
Bð1  mÞ 
T abc ¼ 2r;n ½ð1  mÞdab r;c þ mðdca r;b þ dcb r;a Þ  4ð1 þ mÞr;a r;b r;c 
4pr2

þ2mðna r;b r;c þ nb r;a r;c Þ þ ð1  mÞð2nc r;a r;b þ nb dbc Þ  ð1  3mÞnc dab ðB:6Þ
1 
eabch ¼ 2m½dha r;b r;c þ dhb r;a r;c þ dca r;b r;h þ dcb r;a r;h 
4pr2

þ ð1  mÞ½dbc dah þ dac dbh  dab dch þ 2dab r;c r;h þð1 þ mÞ½2dch r;a r;b  8r;a r;b r;c r;h  ðB:7Þ
 Supriyono, M.H. Aliabadi / International Journal of Solids and Structures 44 (2007) 1038–1059 1055

Appendix C

C.1. Particular solutions for plate bending (Wen et al., 2005)

Governing equation for shear deformable plate bending problem can be written as
^ ¼ Heu
w ðC:1Þ
T T
^ s ¼ f^
where particular solutions of displacement w ^ 3 g ; e ¼ fe1 ; e2 ; e3 g is arbitrary constant vector and
^ 2; w
w1 ; w
the body force in Eq. (25) can be written as
qi ðrÞ ¼ f ðrÞei ðC:2Þ
The components of matrix H are
o2
H ab ¼ 2dab r4  ½ð1 þ mÞr2 þ ð1  mÞk2 
oxa oxb
o
H 3a ¼ H a3 ¼ ð1  mÞðr2  k2 Þ
oxa
and
H 33 ¼ ðr2  k2 Þ½2r2  ð1  mÞk2 =k2 ðC:3Þ
The function u can be defined from Eq. (C.2) such that
Dð1  mÞðr2  k2 Þr4 u þ F ðrÞ ¼ 0 ðC:4Þ
If e1 = 0, e2 = 0 and e3 = 1, the particular solution for plate bending can be written as
1 ou
^ ma ¼ 
w
D oxa
and
1
^ m3 ¼
w ½2r2 u  ð1  mÞk2 w ðC:5Þ
ð1  mÞDk2
where
r4 uðrÞ þ F ðrÞ ¼ 0 ðC:6Þ
The particular solutions of moment and shear force can be determined from shear deformable plate bending
stress resultant-displacement relationship. The tractions on the boundary can be obtained by
b ab nb ;
^pma ¼ M b a na
pm3 ¼ Q
^ ðC:7Þ
If radial basis function F(r) = 1 + r, The function u(r) can be solved from Eq. (C.6)
 4 
r r5
uðrÞ ¼  þ ðC:8Þ
64 225
and the rotations and deflection can be deduced
 
m 1 r x1 r 2
^ 31 ¼ 
w þ
16 45 D
  ðC:9Þ
m 1 r x2 r 2
^ 32 ¼ 
w þ
16 45 D
and
   
1 2r r2 1 r r4
^ m33
w ¼ þ þ þ
2 9 ð1  mÞk2 D 64 225 D
1056 Supriyono, M.H. Aliabadi / International Journal of Solids and Structures 44 (2007) 1038–1059

The particular solutions of moments M b b can be determined from shear deformable

b ab and shear forces Q
plate bending stress resultant–displacement relationship to give
   2 
b m 1 r 2 2 r r3
M 311 ¼  þ ðx1 þ mx2 Þ þ ð1 þ mÞ þ
8 15 16 45
 
Mb m ¼ ð1 þ mÞ 1 þ r ðx1 x2 Þ
312
8 15
   2  ðC:10Þ
b m 1 r 2 2 r r3
M 322 ¼  þ ðmx1 þ x2 Þ þ ð1 þ mÞ þ
8 15 16 45
 
b m ¼  x1 1 þ 2r
Q 31
2 3
and
 
b m x2 2r
Q 32 ¼  1þ
2 3
and the tractions on the boundary can be obtained from relationships in Eq. (C.8).
For the derivative of function F,a = xa/r, the solution ua(r) can be found
r 3 xa
ua ðrÞ ¼  ðC:11Þ
45
and particular solutions w ^ amk are
r
^ m11 ¼ ð3x21 þ r2 Þ
w
45D ðC:12Þ
x1 x2 r
^ m12 ¼ 
w
15D
and
rx1
^ 13 ¼ ½30  ð1  mÞk2 r2 
w
45ð1  mÞk2 D
and the particular solutions of moments M b b are
b ab and shear forces Q
 2   2 
b m ¼  x1
M
x1 x
þ 3r þ m 2 þ r
111
15 r r
 2 
x
b m ¼ ð1  mÞ 2 1 þ r x
M 112
15 r
  2   2  ðC:13Þ
b m x 1 x 1 x2
M 122 ¼  m þ 3r þ þr
15 r r
 2 
b m ¼  1 x1 þ r
Q 11
3 r
and
b m ¼  1 x1 x2
Q 12
3 r
for a = 1, and
x1 x2 r
^ m21 ¼ 
w
15D ðC:14Þ
r
^ m21 ¼ ð3x22 þ r2 Þ
w
45D
and
rx2
^ m23 ¼ ½30  ð1  mÞk2 r2 
w
45ð1  mÞk2 D
 Supriyono, M.H. Aliabadi / International Journal of Solids and Structures 44 (2007) 1038–1059 1057

and the particular solutions of moments M b b are

b ab and shear forces Q
 2   2 
b m ¼  x2
M
x1 x
þ r þ m 2 þ 3r
211
15 r r
 2 
x
b m ¼ ð1  mÞ 1 2 þ r x
M 212
15 r
  2   2  ðC:15Þ
b m ¼  2 m 1 þ r þ x2 þ 3r
M
x x
222
15 r r
1
bm ¼  1 2 x x
Q 21
3 r
and
 2 
b m ¼  1 x2 þ r
Q 22
3 r
for a = 2.

C.2. Particular solutions for 2D plane stress (Wen et al., 2005)

An expression for displacement particular solution ^umca can be obtained in polar coordinates with the use of
the Galerkin vector Gab as
1þm c
umca ðrÞ ¼ Gcba;cc ðrÞ 
^ Gca;bc ðrÞ ðC:16Þ
2
where Gab satisfies
2 xc
r4 Gcba þ dcb ¼ 0 ðC:17Þ
ð1  mÞB r

and a solution is determined by

r 3 xc
Gcba ¼  da;b ðC:18Þ
45ð1  vÞB

Substituting Eq. (C.17) into Eq. (C.15), then the particular solution for displacements can be arranged as
  
2 rx1 1 þ m x31
um11 ¼ 
^  þ 3x1 r
ð1  mÞB 3 30 r

and
 2 
ð1 þ mÞ x1 x2
um12
^ ¼ þ x2 r ðC:19Þ
15ð1  mÞB r
and using strain–displacement relationships for 2D plane stress, the strains are obtained as
 2   4 
m 2 x1 r 1þm x1 6x21
^e111 ¼  þ   3þ þ 3r
ð1  mÞ r 3 30 r r
  3 
2 x1 x2 1 þ m x x2 3x1 x2
^em112 ¼    13 þ
ð1  mÞ 6r 30 r r
and
 2 2 
2 1þm x1 x2
^em122 ¼  3 þ 2r ðC:20Þ
ð1  mÞ 30 r
1058 Supriyono, M.H. Aliabadi / International Journal of Solids and Structures 44 (2007) 1038–1059

The particular solution for membrane stress resultant can be derived by substituting Eq. (C.19) into the
stress resultant–strain relationships for 2D plane stress to give
^ m ¼ B½ð1  mÞ^e1 þ m^e1 
N 111 m11 maa
^ m ¼ Bð1  mÞ^e1
N 112 m12

and
^ m ¼ B½ð1  mÞ^e1 þ m^e1 
N ðC:21Þ
122 m22 maa

and the particular solutions of traction are obtained from

^ m nb
^tm ¼ N ðC:22Þ
1a 1ab

In the same way, displacement particular solutions u^m2a can be obtained as follows:
 2 
m ð1 þ mÞ x2 x1
u21 ¼
^ þ x1 r
15ð1  mÞB r
and
  
2 rx2 1 þ m x32
um22 ¼ 
^  þ 3x2 r ðC:23Þ
ð1  mÞB 3 30 r
and the strains are
 2 2 
2 1þm xx
^em211 ¼  1 3 2 þ 2r
ð1  mÞ 30 r
  3 
2 x1 x2 1 þ m x2 x1 3x1 x2
^em212 ¼   3 þ
ð1  mÞ 6r 30 r r
and
   4 
2 x22 r 1þm x2 6x22
^em222 ¼  þ   3þ þ 3r ðC:24Þ
ð1  mÞ r 3 30 r r
The particular solution for membrane stress resultant are
 
N^ m ¼ B ð1  mÞ^e2 þ m^e2
211 m11 maa

N^ m ¼ Bð1  mÞ^e2
212 m12

and
 
^ m ¼ B ð1  mÞ^e2 þ m^e2
N ðC:25Þ
222 m22 maa

and finally the particular solutions of traction are obtained from

^ m nb
^tm2a ¼ N ðC:26Þ
2ab

With approximation function f(r) = 1 + r, the particular solutions of displacement and traction are
  
1  2v0 2 1 10v0 3
umcb ¼ 
^ r ;b r ;c r  3  d bc  r; r;
b c r ðC:27Þ
½5  4v0 G 30½1  v0 G 3
1  2v0 1
umcb;h ¼ 
^ 0
ðdch r;b þ dbh r;c Þr  ½ð9  10v0 Þdbc r;h  dch r;b  dbh r;c  r;b r;c r;h r2 ðC:28Þ
½5  4v G 30½1  v0 G
 
2ð1  2v0 Þ 1 þ v0 1 1 or
^tmcb ¼ r;c nb þ r;b nc þ dcb r
5  4v 0 1  2v 0 2 2 on
 
1 0 0 0 or 2
 ð4  5v Þr;b nc  ð1  5v Þr;c nb þ fð4  5v Þdcb  r;c r;b g r ðC:29Þ
15ð1  v0 Þ on
where G ¼ ð1  mÞB=2; v0 ¼ 1þm
m
.
 Supriyono, M.H. Aliabadi / International Journal of Solids and Structures 44 (2007) 1038–1059 1059

References

Aliabadi, M.H., 2002. The Boundary Element Method, Applications to Solids and Structures. Wiley, Chichester.
Baiz, P.M., Aliabadi, M.H., 2005. Large deflection analysis of shear deformable shallow shells by the field boundary element method. In:
6th Int. Conf. on Boundary Element Tech.
Brebbia, C.A., Nardini, D., 1983. Dynamics analysis in solid mechanics by alternative boundary element procedure. Int. J. Soil Dyn.
Earthquake Eng. 2, 228.
Kamiya, N., Sawaki, Y., 1988. The plate bending analysis by the dual reciprocity boundary elements. Eng. Anal. 5, 46.
Karam, V.J., Telles, J.C.F., 1998. Nonlinear material analysis of Reissner’s plates. In: Plate Bending Analysis with Boundary Element.
Computational Mechanics Publications, Southampton, pp. 127–163.
Lei, X.Y., Huang, M.K., Wang, X.X., 1990. Geometrically nonlinear analysis of a Reissner’s type plate by boundary element method.
Comput. Struct. 37 (6), 911–916.
Naghdi, P.M., 1956. On the theory of thin elastic shells. Quart. Appl. Math. 14, 369–380.
Purbolaksono, J., Aliabadi, M.H., 2005. Large deformation of shear deformable plate by boundary element method. J. Eng. Math. 51,
211–230.
Reissner, E., 1950. On a variational theorem in elasticity. J. Math. Phys. 29, 90–95.
Ribeiro, G.O., Venturini, W.S., 1998. Elastoplastic analysis of Reissner’s plate using the boundary element method. In: Plate Bending
Analysis with Boundary Element. Computational Mechanics Publications, Southampton, pp. 101–125.
Sawaki, Y., Takeuchi, K., Kamiya, N., 1990. Nonlinear bending analysis without domain-cell discretization. Eng. Anal. 7, 130.
Silva, N.A., Venturini, W.S., 1985. Dual reciprocity process applied to solve bending plate on elastic foundation. In: Brebbia, C.A. (Ed.),
Boundary Element X, Stress Analysis. Springer-Verlag, New York, p. 95.
Stricklin, J.A., Haisler, W.E., von Riesemann, W.A., 1972. Computation and solution procedure for nonlinear analysis by combined finite
element-finite difference methods. Comput. Struct. 2, 955–974.
Supriyono, Aliabadi, M.H., 2006. Boundary element method for shear deformable plates with combined geometric and material
nonlinearities. Eng. Anal. Bound. Elem. 30, 31–42.
Wang, W., Ji, X., Tanka, M.A., 1950. A dual reciprocity boundary element approach for problems of large deflection of thin elastic plates.
Comput. Mech. 26, 58–65.
Wen, P.H., Aliabadi, M.H., Young, A., 2005. Large deflection analysis of Reissner plate by boundary element method. Comput. Struct.
83, 870–879.