Computational Mechanics 28 (2002) 275281
Springer-Verlag 2002

DOI 10.1007/s00466-001-0287-6

Stiffened plate bending analysis by the boundary element method

G. R. Fernandes, W. S. Venturini

275
Abstract In this work, the plate bending formulation of publications, have pointed out the capability of the
the boundary element method (BEM) based on the Kir- method for modelling plates in bending, mentioning
chhoffs hypothesis, is extended to the analysis of stiffened further accuracy and reliability.
elements usually present in building floor structures. In order to use BEM to analyse more complex plates,
Particular integral representations are derived to take e.g., stiffened plates of building floor structures, one has to
directly into account the interactions between the beams extend the BEM formulations to take into consideration
forming grid and surface elements. Equilibrium and arbitrarily displayed beams, general boundary and internal
compatibility conditions are automatically imposed by the constraints and several kinds of transversal loads acting
integral equations, which treat this composite structure as over the plate surface or part of it. Along these lines, Song
a single body. Two possible procedures are shown for [6], Hartmann and Zotemantel [7] have presented inter-
dealing with plate domain stiffened by beams. In the first, esting approaches, discussing in detail displacement
the beam element is considered as a stiffer region re- restrictions at internal points and the use of hermitian
quiring therefore the discretization of two internal lines interpolations. More recently, Oliveira Neto and Paiva [8]
with two unknowns per node. In the second scheme, the have shown a BEM/FEM for analysing building floor
number of degrees of freedom along the interface is structures.
reduced by two by assuming that the cross-section motion While BEM is strongly recommended for plate bending
is defined by three independent components only. analysis, in which internal force and displacement fields
are always accurately modelled, the natural choice to solve
Keywords Plate bending, Boundary elements, Building
the building floor structures is the BEM/FEM combina-
floor structures
tions. Boundary elements are recommended to deal with
1 plate elements, which are combined together by enforcing
Introduction equilibrium and compatibility conditions along the inter-
The boundary element method (BEM) is already a well- faces. The FEM is used to model beam elements. However,
established numerical technique to deal with an enormous for complex floors, characterized by a large number of
number of complex engineering problems. Among them, connected beams and many plate regions of different
analysis of plate bending problems has proved to offer a thickness, the number of degrees of freedom rapidly
particularly adequate field of applications for that tech- increases and the solution accuracy diminishes.
nique. The BEM is suitable for evaluating internal force To overcome these difficulties, Venturini and Paiva [9]
concentrations due to loads distributed over small regions, have proposed a scheme to deal with zoned domains
which very often occur in plate bending analysis. More- without dividing them into sub-regions. In their formu-
over, BEM can deal with deflections, slopes, moments, and lations, only the displacements along the interfaces require
shear forces, approximating them by using the same order approximation. Equilibrium conditions are automatically
polynomials. Thus, shear forces are much better evaluated satisfied, so no approximation of these values is required
when compared with other numerical methods; they de- along the interfaces. In addition, the number of degrees of
pend only on the adopted boundary value approximation. freedom at interface nodes is divided by two. A given in-
The first works discussing the use of direct boundary tegral representation of displacements can be easily writ-
element formulation, in conjunction with Kirchhoffs ten for the whole body. This technique has been extended
theory, are of Bezine [1], Stern [2] and Tottenhan [3]. It is to potential and 2D elastic problems [10], preserving the
also important to mention some previous studies dealing original characteristics: reduction of degrees of freedom
with plate bending problems in the context of indirect and increasing result accuracy. More recently, following
methods [4, 5]. These, as well as several other more recent the same ideas, Fernandes et al. [11] have presented a
plate-bending BEM formulation to deal with a varying
thickness problem.
Received 6 November 2000 Here, this formulation is modified to consider plate
domains stiffened by beams without combining the
G. R. Fernandes, W. S. Venturini (&)
Sao Carlos School of Engineering, algebraic equations written separately for each problem.
University of Sao Paulo, Av. Trabalhador Sao-carlense, The beams are treated as regions of different thicknesses,
400, 13566590 Sao Carlos, Brazil which can be very narrow, therefore representing the
e-mails: [email protected], [email protected] stiffness variation introduced by this kind of stiffener.
 Their rigidities are automatically taken into consideration The generalized internal force displacement relations,
by the global integral representations of displacements and

rotations, in which the stiffener rigidity influences are
mij
D mdij w;kk 1
mw;ij 5
given by two-line integrals along the beam sides. To obtain qi
Dw;jji 6
accurately the algebraic equations, quasi-singular integral
schemes are required to perform the integral along the The effective shear force,
beam elements leading to stable numerical solutions. The Vn qn omns =os 7
formulation is further modified by assuming simple dis-
placement approximations in the direction perpendicular where (n, s) are the local co-ordinate system, with n and
to the beam axes. Displacements in this direction are as- s referred to the boundary normal and tangential
sumed to be linear and, therefore, given in terms of two directions, respectively; no summation is implied.
276 The problem definition is then completed by assuming
node values: the deflection and the rotation at the beam
skeleton point. These assumptions reduce strongly the the following boundary conditions over C: ui ui on C1
number of degrees of freedom of the entire problem. (generalized displacements, deflections, and rotations) and
Finally, some numerical examples are presented to pi pi on C2 (generalized tractions, normal bending
illustrate the accuracy of the results, as well as the moments, and effective shear forces), where C1 [ C2 C.
problem size reduction in terms of degree of freedom.
3
2 Integral representations
Basic equations As usual, the integral representations of deflections,
Let us consider a flat plate of thickness h, referred to a slopes, and generalized internal forces for Kirchhoffs
Cartesian system of co-ordinates with axes x1 and x2 lying plates can be derived from a reciprocity relation written in
on its middle surface and axis x3 perpendicular to that terms of bending moments and curvatures of two inde-
plane. The plate domain is denoted by X (see Fig. 1), while pendent mechanical states. The first state is represented by
its boundary is represented by C. It is assumed that a the actual plate bending problem valid over the domain X,
distributed load g is acting in the x3 direction on the plate for which curvatures w;ij q and moments mij q at a field
midplane, with no distributed external moments. point q are defined, as well as the associated boundary
For this plate, the following basic relationships are values: two generalized displacements, ui Q, and two
defined: generalized tractions, pi Q, referred to a boundary field
point Q. The second elastic state is obtained from the
Equilibrium equations in terms of internal forces: fundamental solution w s; q, the deflection at the field
mij ;j
qi 0 1 point q due to a unit load applied at a source point s. From
this state, fundamental values for curvatures, w;ij s; q, and
qi ;i g 0 2 moments, mij s; q, can be easily derived (1,2,3,11). This
where mij are bending and twisting moments, while qi elastic state is defined over the infinite space of domain
represents shear forces, with subscripts taken in the X1 that contains X (see Fig. 1).
range i; j f1; 2g. In order to obtain all the integral representations, the
following reciprocity relation is easily found for the
The plate bending differential equation, constant rigidity case:
Z Z
mij;ij g 0 3
mij s; qw;ij qdXq mij qwij s; qdXq
or X X
g 8
w;iijj i; j 1; 2 4
D
The reciprocity relation (8) is valid only for plates exhib-
3 2
where D Eh =1
m is the flexural rigidity and iting constant rigidity D, but its extension to a general case
w;iijj r2 w, is the Laplacian operator. is simple [11, 12]. Let us consider a plate element now
characterized by exhibiting variable thickness, i.e.,
t tq, which gives corresponding rigidity D Dq.
Following the same steps to derive the reciprocity relations
(8), but now assuming D Dq variable (continuous or
abrupt variation), one can easily obtain,
Z
mij qw;ij s; qdXq
X
Z
Dq 
m s; qw;ij qdXq 9
Do ij
X

where Do is a reference value; For simplicity, it will be

Fig. 1. Plate domain assumed as the plate rigidity at the collocation point Ds.
 From Eq. (9), one can derive the deflection represen- Figure 2 shows the case of a zoned domain with a
tation for the general case of varying rigidity plates with or narrow region, which can be degenerated to represent the
without abrupt rigidity changes. Thus, assuming that the behaviour of a beam element. For simplicity, single sub-
rigidity can vary over the domain, one can integrate scripts are adopted to represent the external boundary.
Eq. (9) by parts to obtain the deflection integral repre- In order to write an appropriate number of equations
sentation, as follows, relating boundary and interface values, the usual way is to
XNs Z  differentiate equation (11) to obtain other integral
1
CSwS DQVn S; Q representation, finding the free terms that may appear due
k1
DS to differentiation of strong singular kernels [13]. The slope
Ck
 integral representation is the usual relationship frequently
oDQ  oDQ 
2 Mns S; Q Mn S; Q wQdCQ adopted to complete the necessary number of relations for 277
os on the corresponding algebraic system of equations. For an
XNs Z
1 ow internal point S of sub-region Xi the slope equations is

DQMn S; Q QdCQ given by,
k1
DS on
Ck
X
Ns Z   

XNc Dk oVn oM ow
DC  ows=om w
n dC
RC S; CwC C D om om on
DS k1 i
C1 Ckj
Z  2 
Z  
1 oDq  o Dq  XNc

2 qi S; q mij S; q wqdXq Dc oRC ow o2 w
DS oxi oxi oxj w Vn
Mn dC
X D om C
C1 i
om onom
Z   C
ow Z  
Vn Qw S; Q
Mn Q S; Q dCQ XNc
on owC ow
C RC g dX
Z C1
om om
X
Nc
Xg
RC CwC S; C gqw S; qdXq
C1
Xg
12
10 where m defines the slope direction.
For boundary and interface points the vector m must be
where the rigidity DS at source point S was assumed as properly defined according to the adopted rotation value.
Do , Ck is the sub-region boundary which may be defined to As is well known in the case of constant thickness plate,
make possible abrupt rigidity variation, while Ns is the boundary and interface values have to be approximated
total number of sub-regions. before algebraic relations can be obtained from the
It is worth noting that only deflection and rotation integral representations (11) and (12). Although using only
values are required at points along the interface. The algebraic deflection relations to find the system of equa-
right-hand side integral over C is only performed over tions has already proved to be an appropriate scheme for
the external boundary. Thus, traction values along the analysing plate bending problems [11], in the present
interface were eliminated, automatically satisfying formulation it is convenient to also use slope representa-
equilibrium conditions. tions for collocations defined along the interfaces. More-
Equation (10) can be particularized for the simple case over, at corners and other boundary and interface nodes
of zoned domain, in which the rigidities are constant over characterized by the presence of a discontinuity, double
the sub-regions, therefore the first and second derivatives nodes are adopted. In this case, an extra equation must be
of the sub-region rigidity in Eq. (10) are assumed null. written to make corner reaction and the corresponding
This particular case has already been discussed in previous deflection independent values. Thus, at those points where
publications where alternative strategies have been used to discontinuity is assumed, five independent equations are
obtain the final integral representation [11, 12]. required: three deflection representations and two slope
Thus, for a collocation point Si belonging to sub-region representations. For nodes along the external boundary,
i and considering only the case of abrupt rigidity changes, two deflection equations are adopted, one for the bound-
Eq. (10) is reduced to:
XNs Z   XNc
Dk   ow Dc 
CwSi Vn w
Mn dC R w
D
k1 i
on D C C
C1 i
Ckj
Z   XNc Z
 ow 
Vn w
Mn dC RC wC gw dX
on C1
C Xg

11
where Di is the rigidity of the sub-region Xi taken as the
reference value and Ckj is the interface between the sub-
region Xk and its adjacent sub-region Xj . Fig. 2. Zoned domain with a narrow region representing a beam
 ary node and the other for an external collocation very Bending and twisting moment integral representations are
close to the boundary. For interface nodes, deflection and obtained from Eq. (13) by simply applying the definition
slope representations are adopted. given in Eq. (5). To complete the internal force values at
Free term definition is another important aspect internal points, one can differentiate Eq. (13) once more
requiring discussion. For simplicity, slope equations are and apply the definition of shear forces of Eq. (6) to find
only written for smooth boundary and interface nodes, its integral representation.
which lead to simple values. However, it is difficult to The derivatives appearing in Eq. (13) and in the shear
avoid computing the free terms of the deflection equation representation must be performed carefully; free terms
for collocations defined at corners. In Fig. 3, the usual always arise when performing derivatives of strong and
cases of external and internal corners for which one needs hypersingular integrals. Relevant information about the
278 to derive the free terms are depicted, while the corre- correct procedure for extending Eq. (13) can be found
sponding values are given in Table 1. elsewhere [13].
In Table 1, Da is the rigidity of the adjacent sub-region,
while b2 and c are internal angles of the corresponding 4
sub-regions, values usual in BEM formulations. Algebraic equations
To complete the necessary integral relations for the As usual for any BEM formulation, the integral
analysis of plate bending problems considering beam representations (11), (12), (13) and the ones written to
inclusion, one can differentiate Eq. (11) twice to obtain compute internal forces can be transformed into algebraic
the curvature integral representation at internal points, as expressions after discretizing the boundary into elements.
follows, In the present case, plate boundary and interface have
Z   been discretized into geometrically linear elements over
o2 wqm X Ns
Dk o2 Vn o2 Mn ow
w
dC which the boundary values have been approximated by
oxi oxj D
k1 m
oxi oxj oxi oxj on quadratic shape functions. The approximations of the
Ck
  boundary and interface values enable us to write algebraic
XNc
Dc o2 Rc representations of deflections, rotations, moments, and
w shear forces for any collocation point taken inside the
c1 m
D oxi oxj c
Z     domain, along the boundary, along the interfaces, and
o2 w o2 ow outside the plate. After selecting an appropriate number of
Vn
Mnn dC algebraic equations, one can assemble a convenient set of
oxi oxj oxi oxj on
C equations to solve the problem in terms of boundary and
X
Nc Z  2  interface values. This set of algebraic equations has been
o2 wc ow
Rc g dX 13 defined by writing two relations per node plus an extra
i1
oxi oxj oxi oxj relation at each internal or external corner. Only the
Xg
deflection representations are defined for boundary nodes
while, for interface boundary nodes, slope representations
are also adopted. The collocations related with boundary
nodes are taken either along the boundary or outside the
body, while for interface nodes the collocations are always
defined along the interface line.
The outside collocations have been placed very close to
the boundary. To guarantee that the algebraic re-
presentations are accurate, a sub-element scheme has been
used to compute the boundary and interface element in-
tegrals. In this case, sub-element lengths must be smaller
than half of the minimum distance of the collocation to
any sub-element point. This scheme has proved to be very
Fig. 3. Internal and external corner type cases precise and guarantees the quality of the final results.
After performing the element integrals, the algebraic set
Table 1. Free terms K(Q) for usual corners in zoned domains of equations reads:
K(Q) Type of corner HU GP T 14
0.5 Smooth external boundary collocation where U contains deflection and rotation boundary and
  node Q interface nodal values, while P contains boundary node
1 Da
2 1 Dj Smooth interface collocation node tractions only; T is the independent vector due to the
Q; coefficient referred to region Wj
  applied loads.
b3 Da Type C3 corner collocation node Q
1 2p Dj
1 This scheme is simple and can be applied to analyse
b2
Da c Type C2 corner collocation node Q; plates stiffened by beam elements accurately. One needs
2p Dj 2p
b2 is referred to the plate with rigidity Dj; only to model the beams as narrow plate regions, which
c is referred to the plate with rigidity Da leads to a system of equations with a reduced number of
b
Corner collocation along external
2p unknowns. There are two important aspects that deserve
boundary; node type C1
special attention when using this scheme. Quasi-singular
 279
Fig. 4. Plate domain stiffened by beams represented by narrow
regions

integrals are necessarily present due to the narrowness of

some regions. They must be performed accurately by using Fig. 5. Assumed deflections and rotations over the beam
either analytical expressions or numerical schemes with an cross-section
appropriate sub-elementation procedure. For practical
building floor applications, performing numerically the the whole displacement filed over these regions. Several
integrals has proved to be accurate enough. other alternatives have been tested, for which more de-
When the beam is too narrow, the algebraic system of grees of freedom to model the intersection zones were
equations becomes unstable: displacement equations of required. It is important to observe that the skeleton line
any two close collocation points are so similar that they nodes and the geometry of the beam cross-sections are the
lead to inaccurate numerical solutions. In this case, to only data required to define the narrow region, i.e., the
improve the quality of the results one has to eliminate all beam. The two close actual interfaces can then be properly
strongly dependent degrees of freedom, eliminating as well generated to perform the integrals. The final system of
the corresponding algebraic relations. The simplest way to equations is now given in terms of four boundary values
reduce the number of unknowns is given by assuming (w; ow=on, Mn and Vn ) and two interfaces values
rigid rotation of the beam cross-section, making deflec- (w; ow=on).
tions and rotations of two opposite points no longer No limit to reduce the integrals to one single line has
independent values. Instead, they are written in terms of been set. The only approximation assumed is the rigid
the values of nodes defined along the beam skeleton lines. rotation given by Eq. (15). This scheme leads to very small
Rotation will be assumed constant over the section, while system of algebraic equations without introducing too
deflections are given by the skeleton values corrected many variable approximations.
according to the distance from this line. Figure 4 defines
narrow regions as representing beams. The cross-section 5
rigid rotation is illustrated in Fig. 5, from which one may Examples
write a simple relation between rotation, deflection values The stiffened plate analysed here and the two discretiza-
taken on the right (Crib ), left (Clib ) interfaces, and along the tions following the schemes described in the previous
skeleton. sections are defined in Fig. 6. One central beam plus two
Thus, according to these kinematic assumptions, de- others placed along the boundary are adopted to stiffen
flections and rotations along both sides (Crib and Clib ) of the this rectangular plate. The two stiffened sides are free,
narrow region of an internal beam read: while the two others are simply supported. The load is
al
ar ac 15a given by the distributed moment Mn 10 kNm/m applied
along the simply supported sides. Plate and beam mate-
wr wc ac b=2; w wc
ac b=2 15b rials are assumed elastic with Young modulus
where wr , w and wc are deflections on the right, left, and E 3 000 000 kN/cm2 and Poisson ration m 0:16. The
along the skeleton, respectively, while ar , a and ac are the plate thickness is hp 8 cm, while the beam heights are
rotations. hb 20 cm. The discretizations are specified by: (a) using
In the case of an external beam, the deflections along the scheme where interface collocations are considered,
the sides (Creb and Cleb ) are also given by the expression leading to 38 elements and 92 nodes (see Fig. 6b) and; (b)
(15b); however, for the rotations one has: defining collocations only along the beam skeleton lines,
leading to 22 elements and 36 nodes (see Fig. 6a). Figures 7
al ar ac 15c
and 8 give the displacements and the moments in x and y
With this approximation only two algebraic relations are directions, respectively. The values were obtained along
required at each beam node. Thus, the collocation point is the central beam (V2 ) by using the two discussed dis-
now placed along the beam skeleton, while wc and ac are cretization schemes. As one can see, the results given in
the assumed unknown values. The intersection regions terms of deflections obtained by using both schemes
between two or more beams were also assumed rigid and, compare very well (Fig. 7). Comparing the moment dis-
consequently, only three displacements are left to describe tribution along the central beam obtained for both cases,
280

Fig. 8. Moments Mx and My along the central beam (V2 )

other two sides are free. The problem has also been solved
by the two proposed schemes. Using the first scheme,
characterized by defining collocations along the interfaces,
275 nodes were defined along 116 boundary and internal
elements. When the second scheme is adopted, the total
Fig. 6. Stiffened plate
number of nodes is reduced to 123 with 53 elements (see
Fig. 9). The results of both analyses in terms of deflections
are displayed in Fig. 10, while Fig. 11 gives the corre-
sponding distribution of moments, Mx and My , computed
along the central beam (V5 ). Again, the computed deflec-
tions using both schemes compare very well. The numer-
ical results in terms of moment components appear to be
quite reasonable for this case. The only difficult appears to
be the accurate computation of the moment values at the
beam intersection. The differences observed at this point
are introduced by the assumed rigid rotations.

6
Conclusions
The BEM formulation for analysing zoned plate-bending
Fig. 7. Deflections along the central beam (V2 ) problem has been extended to deal with plate stiffened
by beams. Beam rigidity is taken into account by
assuming narrow sub-regions, without dividing the
very small differences are observed (Fig. 8). Thus, the rigid stiffened plate into beam and plate elements. Equilibrium
rotation hypothesis assumed for the second scheme in- and compatibility conditions are automatically guaran-
troduces negligible differences in the actual beam stiffness teed by the global integral equations. Two schemes
as well as in the whole structure. Assuming that the results regarding the displacement approximation assumed
are enough adequate for analysing stiffened building floors along the beams were proposed. The formulation avoids
for instance, the second scheme is recommended since it unnecessary approximations usually present when treat-
requires a reduced number of algebraic equations, conse- ing this problem with the standard sub-region technique,
quently reducing the computations. which reduces the number of unknowns and increases
To run the second example, we have borrowed several the accuracy of the results. The rigid cross-section
data from the previously analysed problem: the applied assumption made for the second scheme reduces even
boundary moment Mn , the Young modulus E, the Pois- further the number of unknowns and preserves the
sons ratio m, the plate thickness hp , and the beam height quality of the results when analysing stiffened plates
hb . Here, the plate is stiffened by beams placed along all often adopted for building floor structures showing their
sides and by two internal ones, as seen in Fig. 9. Boundary adequacy for dealing with practical problems, thus
values given in terms of moments are prescribed along the encouraging the BEM researchers to go ahead with this
simply supported sides connected to beams V1 and V3 ; the kind of formulation.
 281

Fig. 9. Stiffened plate

References
1. Bezine GP (1978) Boundary integral formulation for plate
flexure with arbitrary boundary conditions. Mech. Res.
Comm. 5(4): 197206
2. Stern MA (1979) A general boundary integral formulation for
the numerical solution of plate bending problems. Int.
J. Solids Struct. 15: 769782
3. Tottenhan H (1979) The boundary element method for plates
and shells. In: Banerjee PK, Butterfield R (eds) Developments
in Boundary Element Methods, pp. 173205
4. Jaswon MA, Maiti M, Symm GT (1967) Numerical bihar-
monic analysis and some applications. Int. J. Solids Struct. 3:
309332
5. Jaswon MA, Maiti M (1968) An integral formulation of plate
bending problems. J. Eng. Math. 2: 8393
6. Song GS, Mukherjee S (1986) Boundary element method
Fig. 10. Deflections along the central beam (V5 ) analysis of bending of elastic plates of arbitrary shape with
general boundary conditions. Eng. Anal. Boundary Elements
3: 3644
7. Hartmann F, Zotemantel R (1986) The direct boundary ele-
ment method in plate bending. Int. J. Num. Meth. Eng.
23(11): 20492069
8. Oliveira Neto L, Paiva JB (1998) A three nodal parameters
BEM formulation for building floor slabs analysis. In: World
Congress on Computational Mechanics, 4., Buenos Aires,
Argentina, Abstracts. Santa Fe, IACM, v. 1, p. 326
9. Venturini WS, Paiva JB (1993) Boundary elements for plate
bending analysis. Eng. Anal. 11: 18
10. Venturini WS (1992) Alternative formulations of the
boundary element method for potential and elastic zoned
problems. Eng. Anal. 9: 203207
11. Chaves EWV, Fernandes GR, Venturini WS (1999) Plate
bending boundary element formulation considering variable
thickness. Eng. Anal. Boundary Elements 23: 405418
12. Paiva JB, Aliabadi MH (1998) Boundary element analysis of
zoned bending plates. Proceedings of IABEM98, Ecole Poly-
technique, Palaiseau (Paris)
13. Frangi A, Guiggiani M (1999) Boundary element analysis of
Kirchhoff plates with direct evaluation of hypersingular
integrals. Int. J. Num. Meth. Eng. 46: 18451863
Fig. 11. Moments Mx and My along the central beam (V5 )