# Nonlinear Dynamic Analysis of Plates Stiffened by Parallel Beams with Deformable Connection

Authors: J. A. Dourakopoulos; E. J. Sapountzakis
Journal: Advances in Civil Engineering
 (2014)

---

## Abstract

In this paper a general solution to the geometrically nonlinear dynamic analysis of plates stiffened by arbitrarily placed parallel beams of arbitrary doubly symmetric cross-section, subjected to dynamic loading, is presented. The plate-beam structure is assumed to undergo moderate large deflections and the nonlinear analysis is carried out by retaining nonlinear terms in the kinematical relations. According to the proposed model, the arbitrarily placed parallel stiffening beams are isolated from the plate by sections in the lower outer surface of the plate, making the hypothesis that the plate and the beams can slip in all directions of the connection without separation and taking into account the arising tractions in all directions at the fictitious interfaces. These tractions are integrated with respect to each half of the interface width resulting in two interface lines, along which the loading of the beams and the additional loading of the plate are defined. Six boundary value problems are formulated and solved using the analog equation method (AEM), a BEM-based method. Both free and forced transverse vibrations are considered and numerical examples with great practical interest are presented demonstrating the effectiveness, wherever possible, the accuracy, and the range of applications of the proposed method.


---

## Body

## 1. Introduction

Stiffened plate panels are structural elements of practical importance in applications such as ship superstructures, bridge decks, and aircraft structures. Stiffening of the plate provides the benefit of added load-carrying capability with a relatively small additional weight penalty, while buckling is prevented especially in case of inplane loading. The unique peculiarities of the aforementioned structures are obtained due to the behavior of the bond between the plate and the beams; however, this bond is the usual reason why these structures are prone to failure. Moreover, since these stiffened plates are frequently subjected to dynamic loading such as air blasts or underwater explosions, a clear understanding of the dynamic response requires development of an efficient dynamic analysis capability.When the deflections of the structure are small, a wide range of linear analysis tools, such as modal analysis, can be used and some analytical results are possible. As the deflections become larger, the induced geometric nonlinearities result in effects that are not observed in linear systems, making therefore the determination of an analytical solution extremely difficult and in most cases impossible. Moreover, having in mind the importance of weight saving in engineering structures, the study of nonlinear effects in the analysis (large deflection analysis) of stiffened plates becomes essential. These nonlinearities result from retaining the squares of the slopes in the strain-displacement relations (intermediate nonlinear theory), avoiding in this way the inaccuracies arising from a linear or a linearized second-order analysis.The linear dynamic behavior of stiffened plates has been widely studied employing the Rayleigh-Ritz method [1–4], the transfer matrix method [5], the finite difference method [6–9], the finite element method [10–12], the finite strip method [13], and the boundary element method (BEM) [14–16]. To the authors’ knowledge, a limited amount of technical literature is also available on the nonlinear dynamic analysis of stiffened plate systems. According to this, the integral equation method has been applied for the large deflection analysis of clamped laterally loaded skew plates with stiffener parallel to the skew directions [17]. Moreover, the finite element method has been widely used for the free large-amplitude flexural vibration of stiffened [18] and laminated plates [19, 20] and for forced vibrations under instantaneous loading taking into account both geometric and material nonlinearities [21]. The spline finite strip method has also been applied for the examination of nonlinear transient vibration of stiffened plates [22] and the Galerkin method has been employed for the solution of discrete equations derived after the application of Fourier series for the deflection of the plate [23, 24]. Nevertheless, in all of the aforementioned research efforts restrictions are encountered, arising from the ignorance of the tractions in the fictitious plate-beams interfaces, of the nonuniform distribution of the interface transverse shear force, or of the nonuniform torsional response of the beams.The behavior of composite slab-and-beam structure is affected significantly by the deformability of the shear connection. Therefore, much research effort has been done concerning the linear dynamic analysis of composite beams with deformable connections. Wu et al examined the free vibration of partial interaction composite beams with axial force and proposed an approximate simple expression to predict the fundamental frequency of the partial-interaction composite members with axial force [26]. Xu and Wu investigated the static, dynamic, and buckling behavior of partial-interaction composite members by taking into account the influences of rotary inertia and shear deformations [27]. Moreover, they also employed the state-space method in order to develop an exact two-dimensional plane stress model of composite beams with interlayer slips [28, 29].On the other hand, less research effort has been made on the forced vibration of partial-interaction composite beams. Girhammar and Pan [30] investigated the dynamic behavior of partial-interaction composite beams without axial force and presented the orthogonality relation of vibration modes. They also obtained the analytical expressions for the transient response of composite beams under impulsive and step loads. Adam et al. [31] decomposed the dynamic response into two parts, the quasi-static component solved by integrating the influence function and the complementary component solved by the mode superposition. Recently, Girhammar et al. [32] employed Hamilton’s principle to derive the governing equations for dynamic problems of composite beams under generalized boundary conditions. At the same time, they studied the forced vibration of composite beams by virtue of mode superposition. However, in the literature, large deflection analysis of partial interaction has received limited attention and focused especially on the static behavior of the structure [33–38].In this paper, a general solution to the geometrically nonlinear dynamic analysis of plates stiffened by arbitrarily placed parallel beams of arbitrary doubly symmetric cross-section with deformable connections subjected to arbitrary dynamic loading is presented. The solution is based on the structural model proposed by Sapountzakis and Mokos in [39] and Dourakopoulos and Sapountzakis in [40], according to which the stiffening beams are isolated from the plate by sections in the lower outer surface of the plate, making the hypothesis that the plate and the beams can slip in all directions of the connection without separation (i.e., uplift is neglected) and taking into account the arising tractions in all directions at the fictitious interfaces. Due to the slip, relative inplane displacements between beam and plate arise, assuming however that they are continuously in contact and taking into account the arising tractions in all directions at the fictitious interfaces. These tractions are integrated with respect to each half of the interface width resulting in two interface lines, along which the loading of the beams and the additional loading of the plate are defined. The unknown distribution of the aforementioned integrated tractions is established by applying continuity conditions in all directions at the two interface lines taking into account their relationship with the interface slip through the shear connector stiffness. Any distribution of connectors in each direction of the interfaces can be handled. The utilization of two interface lines for each beam enables the nonuniform torsional response of the beams to be taken into account as the angle of twist is indirectly equated with the corresponding plate slope. Six boundary value problems are formulated and solved using the analog equation method (AEM) [41], a BEM-based method. The essential features and novel aspects of the present formulation are summarized as follows.(i)
The adopted model permits the evaluation of the inplane shear forces (longitudinal and transverse) at the interfaces in both directions, taking into account the influence of interface slip, the knowledge of which is very important in the design of shear connectors in plate-beam structures.
(ii)
Utilization of two interface lines permits the nonuniform distribution of the inplane forces along the interface width to be taken into account.
(iii)
Both free and forced transverse vibrations are considered taking into account geometric nonlinearities.
(iv)
The stiffened plate is of arbitrary shape and is subjected to arbitrary dynamic loading, while both the number and the placement of the parallel beams are also arbitrary (eccentric beams are included).
(v)
The cross-section of the stiffening beams is an arbitrary doubly symmetric thin or thick-walled one. The formulation does not stand on the assumption of a thin-walled structure according to which torsional and warping rigidities can be evaluated employing closed-form expressions without significant error as long as the cross-section thickness is small [42]. Therefore, in this paper, the cross-section’s torsional and warping rigidities are evaluated “exactly” in a numerical sense employing the BEM (requiring only boundary discretization for the cross-sectional analysis).
(vi)
The plate and the beams are supported by the most general boundary conditions including elastic support or restraint.
(vii)
The nonuniform torsion in which the stiffening beams are subjected is taken into account by solving the corresponding problem and by comprehending the arising twisting and warping in the corresponding displacement continuity conditions. The distributed warping moment arising from the nonuniform distribution of longitudinal inplane forces is also taken into account.
(viii)
Contrary to previous research efforts where the numerical analysis is based on BEM using a lumped mass assumption model after evaluating the flexibility matrix at the mass nodal points, in this work, a distributed mass model is employed.Based on the numerical solution developed, several examples with great practical interest are presented, demonstrating the effectiveness, wherever possible, the accuracy, and the range of applications of the proposed method.


## 2. Statement of the Problem

Let us consider a thin plate of homogeneous, isotropic, and linearly elastic material with modulus of elasticityE
p, volume mass density ρ
p, shear modulus G
p, and Poisson’s ratio ν
p having constant thickness h
p and occupying the two-dimensional multiple connected region Ω of the x
y plane bounded by the piecewise Γ
j
(
j
=
0,1
,
2
,
…
,
K
) boundary curves. The plate is stiffened by a set of i
=
1,2
,
…
,
I arbitrarily placed parallel beams of arbitrary doubly symmetric cross-section of area A
b
i and length l
i. The material of the beams is considered to be homogeneous, isotropic, and linearly elastic with modulus of elasticity E
b
i, mass density ρ
b
i, shear modulus G
b
i, and Poisson’s ratio ν
b
i. Therefore, material nonlinearities are not considered in the following analysis.For the sake of convenience thex axis is taken parallel to the beams of length l
i which may have either internal or boundary point supports. The stiffened plate is subjected to the arbitrary lateral dynamic load g
=
g
(
x
,
t
), x
:
{
x
,
y
},  t
≥
0. Owing to this loading, the plate and the beams can slip in all directions of the connection without separation (i.e., uplift is neglected). For the analysis of the aforementioned problem, a global coordinate system O
x
y for the analysis of the plate and local coordinate ones C
i
x
i
y
i corresponding to the centroid axes of each beam are employed (Figure 1).Figure 1
Two-dimensional regionΩ occupied by the plate.The solution of the problem at hand is approached employing the improved model proposed by Sapountzakis and Mokos in [39] and Dourakopoulos and Sapountzakis in [40]. According to this model, the stiffening beams are isolated from the plate by sections in its lower outer surface, taking into account the arising tractions at the fictitious interfaces (Figure 2). Integration of these tractions along each half of the width of the ith beam results in line forces per unit length in all directions in two interface lines, which are denoted by q
x
j
i, q
y
j
i, and q
z
j
i (j
=
1,2), encountering in this way the nonuniform distribution of the interface transverse shear forces q
y
i, which in previous models [43] was ignored. The aforementioned integrated tractions result in the loading of the ith beam and the additional loading of the plate. Their distribution is unknown and can be established by imposing displacement continuity conditions in all directions along the two interface lines, enabling in this way the nonuniform torsional response of the beams to be taken into account. Thus, the arising additional loading (Figure 3) at the middle surface of the plate and the loading along the centroid axis (coinciding with the shear center axis), of each beam, can be summarized as follows.Thin elastic plate stiffened by beams (a) and isolation of the beams from the plate (b).
(a)
(b)Figure 3
Structural model and directions of the additional loading of the plate and theith beam.(a) In the Plate (at the Traces of the Two Interface Lines
j
=
1,2
of the
i
th Plate-Beam Interface)
(i)
An inplane line body forceq
x
j
i at the middle surface of the plate.
(ii)
An inplane line body forceq
y
j
i at the middle surface of the plate.
(iii)
A lateral line loadq
z
j
i.
(iv)
A lateral line load∂
m
p
y
j
i
/
∂
x due to the eccentricity of the component q
x
j
i from the middle surface of the plate. m
p
y
j
i
=
q
x
j
i
h
p
/
2 is the bending moment.
(v)
A lateral line load∂
m
p
x
j
i
/
∂
y due to the eccentricity of the component q
y
j
i from the middle surface of the plate. m
p
x
j
i
=
q
y
j
i
h
p
/
2 is the bending moment.(b) In Each (
i
th) Beam (
C
i
x
i
y
i
z
i
System of Axes)
(i)
An axially distributed line loadq
x
j
i along the beam centroid axis C
i
x
i.
(ii)
A transversely distributed line loadq
y
j
i along the beam centroid axis C
i
x
i.
(iii)
A perpendicularly distributed line loadq
z
j
i along the beam centroid axis C
i
x
i.
(iv)
A distributed bending momentm
b
y
j
i
=
q
x
j
i
e
z
j
i along C
i
y
i local beam centroid axis due to the eccentricities e
z
j
i of the components q
x
j
i from the beam centroid axis. e
z
1
i
=
e
z
2
i
=
-
h
b
i
/
2 are the eccentricities.
(v)
A distributed bending momentm
b
z
j
i
=
-
q
x
j
i
e
y
j
i along C
i
z
i local beam centroid axis due to the eccentricities e
y
j
i of the components q
x
j
i from the beam centroid axis. e
y
1
i
=
-
b
f
i
/
4 and e
y
2
i
=
b
f
i
/
4 are the eccentricities.
(vi)
A distributed twisting momentm
b
x
j
i
=
q
z
j
i
e
y
j
i
-
q
y
j
i
e
z
j
i along C
i
x
i local beam shear center axis due to the eccentricities e
z
j
i and e
y
j
i of the components q
y
j
i and q
z
j
i from the beam shear center axis, respectively. e
z
1
i
=
e
z
2
i
=
-
h
b
i
/
2 and e
y
1
i
=
-
b
f
i
/
4 and e
y
2
i
=
b
f
i
/
4 are the eccentricities.
(vii)
A distributed warping momentm
b
w
j
i
=
-
q
x
j
i
(
φ
S
P
i
)
f
j,along C
i
x
i local beam shear center axis, was ignored in previous models [39, 43]. (
φ
S
P
i
)
f
j is the value of the primary warping function φ
S
P
i with respect to the shear center of the beam cross-section (coinciding with its centroid) at the point of the jth interface line of the ith plate-beam interface.On the basis of the above considerations the response of the plate and the beams may be described by the following boundary value problems.(a) For the Plate.The analysis of the plate is based on the Von Kármán plate theory, according to which the deflection of the plate cannot be regarded as small as compared to the plate thickness, while it remains small in comparison with the rest dimensions of the plate. Due to this assumption, geometrical nonlinearities should be taken into account and the displacement field of an arbitrary point of the plate, as implied by the Kirchhoff hypothesis, is given as(1a)
u
-
p
(
x
,
y
,
z
,
t
)
=
u
p
(
x
,
y
,
t
)
-
z
∂
w
p
(
x
,
y
,
t
)
∂
x
,
(1b)
v
-
p
(
x
,
y
,
z
,
t
)
=
v
p
(
x
,
y
,
t
)
-
z
∂
w
p
(
x
,
y
,
t
)
∂
y
,
(1c)
w
-
p
(
x
,
y
,
z
,
t
)
=
w
p
(
x
,
y
,
t
)
,where u
-
p, v
-
p, w
-
p, x
:
{
x
,
y
}, t
≥
0 are the time dependent inplane and transverse displacement components of an arbitrary point of the plate and u
p
=
u
p
(
x
,
t
), v
p
=
v
p
(
x
,
t
), and w
p
=
w
p
(
x
,
t
),x
:
{
x
,
y
}, t
≥
0 are the corresponding components of a point at its middle surface. Employing the strain-displacement relations of the three-dimensional elasticity for moderate large displacements [44, 45], the strain components can be written as(2a)
ε
x
x
=
∂
u
-
p
∂
x
+
1
2
(
∂
w
-
p
∂
x
)
2
,
(2b)
ε
y
y
=
∂
v
-
p
∂
y
+
1
2
(
∂
w
-
p
∂
y
)
2
,
(2c)
γ
x
y
=
∂
u
-
p
∂
y
+
∂
v
-
p
∂
x
+
∂
w
-
p
∂
x
∂
w
-
p
∂
y
,
(2d)
ε
z
z
=
γ
x
z
=
γ
y
z
=
0
.Substituting (1a)–(1c) and (2a)–(2d) to the stress-strain relations defined by the Hooke’s law
(3)
{
S
x
x
S
y
y
S
x
y
}
=
[
E
p
(
1
-
ν
p
)
2
E
p
ν
p
(
1
-
ν
p
)
2
0
E
p
ν
p
(
1
-
ν
p
)
2
E
p
(
1
-
ν
p
)
2
0
0
0
E
p
2
(
1
+
ν
p
)
]
{
ε
x
x
ε
y
y
γ
x
y
}
the nonvanishing components of the second Piola-Kirchhoff stress tensor are obtained as(4a)
S
x
x
=
E
p
(
1
-
ν
p
2
)
×
{
[
∂
u
p
∂
x
-
z
∂
2
w
p
∂
x
2
+
1
2
(
∂
w
p
∂
x
)
2
]
hhhhhhh
+
ν
p
[
∂
v
p
∂
y
-
z
∂
2
w
p
∂
y
2
+
1
2
(
∂
w
p
∂
y
)
2
]
}
,
(4b)
S
y
y
=
E
p
(
1
-
ν
p
2
)
{
ν
p
[
∂
u
p
∂
x
-
z
∂
2
w
p
∂
x
2
+
1
2
(
∂
w
p
∂
x
)
2
]
hhhhhhhhhhhh
+
[
∂
v
p
∂
y
-
z
∂
2
w
p
∂
y
2
+
1
2
(
∂
w
p
∂
y
)
2
]
}
,
(4c)
S
x
y
=
E
p
2
(
1
+
ν
p
)
(
∂
u
p
∂
y
+
∂
v
p
∂
x
-
2
z
∂
2
w
p
∂
x
∂
y
+
∂
w
p
∂
x
∂
w
p
∂
y
)
.Subsequently, integrating the stress components over the plate thickness, the stress resultants acting on the plate are written as(5a)
N
p
x
=
C
[
∂
u
p
∂
x
+
ν
p
∂
v
p
∂
x
+
1
2
(
∂
w
p
∂
x
)
2
+
1
2
ν
p
(
∂
w
p
∂
y
)
2
]
,
(5b)
N
p
y
=
C
[
∂
v
p
∂
x
+
ν
p
∂
u
p
∂
x
+
1
2
(
∂
w
p
∂
y
)
2
+
1
2
ν
p
(
∂
w
p
∂
x
)
2
]
,
(5c)
N
p
x
y
=
C
1
-
ν
p
2
(
∂
u
p
∂
y
+
∂
v
p
∂
x
+
∂
w
p
∂
y
∂
w
p
∂
x
)
,
(5d)
M
p
x
=
-
D
(
∂
2
w
p
∂
x
2
+
ν
p
∂
2
w
p
∂
y
2
)
,
(5e)
M
p
y
=
-
D
(
∂
2
w
p
∂
y
2
+
ν
p
∂
2
w
p
∂
x
2
)
,
(5f)
M
p
x
y
=
-
D
(
1
-
ν
p
)
∂
2
w
p
∂
x
∂
y
,where C
=
E
p
h
p
/
(
1
-
ν
p
2
) and D
=
E
p
h
p
3
/
12
(
1
-
ν
p
2
) are the membrane and bending rigidities of the plate, respectively.On the basis of Hamilton’s principle, the system of partial differential equations of motion of the plate in terms of the stress resultants is obtained as(6a)
∂
N
p
x
∂
x
+
∂
N
p
x
y
∂
y
-
ρ
p
h
p
u
¨
p
=
∑
i
=
1
I
(
∑
j
=
1
2
q
x
j
i
δ
(
y
-
y
j
)
)
,
(6b)
∂
N
p
y
∂
y
+
∂
N
p
x
y
∂
x
-
ρ
p
h
p
v
¨
p
=
∑
i
=
1
I
(
∑
j
=
1
2
q
y
j
i
δ
(
y
-
y
j
)
)
,
(6c)
-
∂
2
M
p
x
∂
x
2
-
2
∂
2
M
p
x
y
∂
x
∂
y
-
∂
2
M
p
y
∂
y
2
-
N
p
x
∂
2
w
p
∂
x
2
-
2
N
p
x
y
∂
2
w
p
∂
x
∂
y
-
N
p
y
∂
2
w
p
∂
y
2
+
ρ
p
h
p
w
¨
p
-
ρ
p
h
p
u
¨
p
∂
w
p
∂
x
-
ρ
p
h
p
v
¨
p
∂
w
p
∂
y
-
ρ
p
h
p
3
12
∂
2
w
¨
p
∂
x
2
-
ρ
p
h
p
3
12
∂
2
w
¨
p
∂
y
2
=
g
-
∑
i
=
1
I
[
∑
j
=
1
2
(
q
z
j
i
+
∂
m
p
x
j
i
∂
y
+
∂
m
p
y
j
i
∂
x
-
q
x
j
i
∂
w
p
j
i
∂
x
hhhhhhhhhhhhhhhh
-
q
y
j
i
∂
w
p
j
i
∂
y
)
δ
(
y
-
y
j
)
∑
j
2
]
,where δ
(
y
-
y
i
) is the Dirac’s delta function in the y direction. Employing relations (5a)–(5f), the governing differential equations (6a)–(6c) in the domain Ω can be expressed in terms of the displacement components as(7a)
G
p
h
p
[
1
+
ν
p
1
-
ν
p
∂
2
w
p
∂
x
∂
y
∂
w
p
∂
y
∇
2
u
p
+
1
+
ν
p
1
-
ν
p
∂
∂
x
(
∂
u
p
∂
x
+
∂
v
p
∂
y
)
hhhhh
+
(
2
1
-
ν
p
∂
2
w
p
∂
x
2
+
∂
2
w
p
∂
y
2
)
∂
w
p
∂
x
hhhhh
+
1
+
ν
p
1
-
ν
p
∂
2
w
p
∂
x
∂
y
∂
w
p
∂
y
]
-
ρ
p
h
p
u
¨
p
=
∑
i
=
1
I
(
∑
j
=
1
2
q
x
j
i
δ
(
y
-
y
j
)
)
,
(7b)
G
p
h
p
[
1
+
ν
p
1
-
ν
p
∂
2
w
p
∂
x
∂
y
∂
w
p
∂
x
∇
2
v
p
+
1
+
ν
p
1
-
ν
p
∂
∂
y
(
∂
u
p
∂
x
+
∂
v
p
∂
y
)
hhhhh
+
(
2
1
-
ν
p
∂
2
w
p
∂
y
2
+
∂
2
w
p
∂
x
2
)
∂
w
p
∂
y
hhhhh
+
1
+
ν
p
1
-
ν
p
∂
2
w
p
∂
x
∂
y
∂
w
p
∂
x
]
-
ρ
p
h
p
v
¨
p
=
∑
i
=
1
I
(
∑
j
=
1
2
q
y
j
i
δ
(
y
-
y
j
)
)
,
(7c)
D
∇
4
w
p
-
C
×
{
[
(
∂
u
p
∂
x
+
1
2
(
∂
w
p
∂
x
)
2
)
hhhhhhhh
+
ν
p
(
∂
v
p
∂
y
+
1
2
(
∂
w
p
∂
y
)
2
)
]
hhhhhhhh
×
∂
2
w
p
∂
x
2
+
(
1
-
ν
p
)
hhhhhhhh
·
(
∂
u
p
∂
y
+
∂
v
p
∂
x
+
∂
w
p
∂
x
∂
w
p
∂
y
)
∂
2
w
p
∂
x
∂
y
hhhhhhhh
+
[
(
∂
v
p
∂
y
+
1
2
(
∂
w
p
∂
y
)
2
)
hhhhhhhhhhhh
+
ν
p
(
∂
u
p
∂
x
+
1
2
(
∂
w
p
∂
x
)
2
)
]
hhhhhhhh
·
∂
2
w
p
∂
y
2
}
+
ρ
p
h
p
w
¨
p
-
ρ
p
h
p
u
¨
p
∂
w
p
∂
x
-
ρ
p
h
p
v
¨
p
∂
w
p
∂
y
-
ρ
p
h
p
3
12
∂
2
w
¨
p
∂
x
2
-
ρ
p
h
p
3
12
∂
2
w
¨
p
∂
y
2
=
g
-
∑
i
=
1
I
[
∑
j
=
1
2
(
q
z
j
i
+
∂
m
p
x
j
i
∂
y
+
∂
m
p
y
j
i
∂
x
-
q
x
j
i
∂
w
p
j
i
∂
x
hhhhhhhhhhhhhhhhh
-
q
y
j
i
∂
w
p
j
i
∂
y
)
δ
(
y
-
y
j
)
∑
j
=
1
2
(
q
z
j
i
+
∂
m
p
x
j
i
∂
y
+
∂
m
p
y
j
i
∂
x
-
q
x
j
i
∂
w
p
j
i
∂
x
]
.The governing differential equations (7a)–(7c) are also subjected to the pertinent boundary conditions of the problem at hand(8a)
a
p
1
u
p
n
+
a
p
2
N
p
n
=
a
p
3
,
(8b)
β
p
1
u
p
t
+
β
p
2
N
p
t
=
β
p
3
,
(8c)
γ
p
1
w
p
+
γ
p
2
R
p
n
=
γ
p
3
,
(8d)
δ
p
1
∂
w
p
∂
n
+
δ
p
2
M
p
n
=
δ
p
3
,
(8e)
ε
1
k
w
p
+
ε
2
k
∥
T
w
p
∥
k
=
ε
3
k
,
ε
2
k
≠
0and to the initial conditions(9a)
w
p
(
x
,
0
)
=
w
p
0
(
x
)
,
(9b)
w
˙
p
(
x
,
0
)
=
w
-
p
0
(
x
)
,where a
p
l, β
p
l, γ
p
l, and δ
p
l
(
l
=
1,2
,
3
) are functions specified at the boundary Γ; ε
l
k
(
l
=
1,2
,
3
) are functions specified at the k corners of the plate; w
p
0
(
x
), w
-
p
0
(
x
), and x
:
{
x
,
y
} are the initial deflection and velocity of the points of the middle surface of the plate; u
p
n, u
p
t and N
p
n, N
p
t are the boundary membrane displacements and forces in the normal and tangential directions to the boundary, respectively; R
p
n and M
p
n are the effective reaction along the boundary and the bending moment normal to it, respectively, which by employing intrinsic coordinates (i.e., the distance along the outward normal n to the boundary and the arc length s) are written as(10a)
R
p
n
=
-
D
[
∂
∂
n
∇
2
w
p
-
(
ν
p
-
1
)
∂
∂
s
(
∂
2
w
p
∂
s
∂
n
-
κ
∂
w
p
∂
s
)
]
+
N
p
n
∂
w
p
∂
n
+
N
p
t
∂
w
p
∂
s
,
(10b)
M
p
n
=
-
D
[
∇
2
w
p
+
(
ν
p
-
1
)
(
∂
2
w
p
∂
s
2
+
κ
∂
w
p
∂
n
)
]in which κ
(
s
) is the curvature of the boundary. Finally, ∥
T
w
p
∥
k is the discontinuity jump of the twisting moment T
w
p at the corner k of the plate, while T
w
p along the boundary is given by the following relation:
(11)
T
w
p
=
D
(
ν
p
-
1
)
(
∂
2
w
p
∂
s
∂
n
-
κ
∂
w
p
∂
s
)
.
The boundary conditions (8a)–(8d) are the most general boundary conditions for the plate problem including also elastic support, while the corner condition (8e) holds for free or transversely elastically restrained cornersk. It is apparent that all types of the conventional boundary conditions can be derived from these equations by specifying appropriately the functions a
p
l, β
p
l, γ
p
l, and δ
p
l (l
=
1,2
,
3) (e.g., for a clamped edge it is a
p
1
=
β
p
1
=
γ
p
1
=
δ
p
1
=
1, a
p
2
=
a
p
3
=
β
p
2
=
β
p
3
=
γ
p
2
=
γ
p
3
=
δ
p
2
=
δ
p
3
=
0).(b) For Each (
i
th) Beam. Each beam undergoes transverse deflection with respect to z
i and y
i axes and axial deformation and nonuniform angle of twist along x
i axis. Based on the Bernoulli theory, the displacement field of an arbitrary point of a cross-section (taking into account moderate large displacements and considering the angle of rotation of twist to have relatively small values) can be derived with respect to those of its centroid as(12a)
u
-
b
i
(
x
i
,
y
i
,
z
i
,
t
)
=
u
b
i
(
x
i
,
t
)
-
y
i
θ
b
z
i
(
x
i
,
t
)
+
z
i
θ
b
y
i
(
x
i
,
t
)
+
∂
θ
b
x
i
∂
x
i
φ
S
P
i
(
y
i
,
z
i
,
t
)
,
(12b)
v
-
b
i
(
x
i
,
y
i
,
z
i
,
t
)
=
v
b
i
(
x
i
,
t
)
-
z
i
θ
b
x
i
(
x
i
,
t
)
,
(12c)
w
-
b
i
(
x
i
,
y
i
,
z
i
,
t
)
=
w
b
i
(
x
i
,
t
)
+
y
i
θ
b
x
i
(
x
i
,
t
)
,
(12d)
θ
b
y
i
(
x
i
,
t
)
=
-
∂
w
b
i
(
x
i
,
t
)
∂
x
i
,
(12e)
θ
b
z
i
(
x
i
,
t
)
=
∂
v
b
i
(
x
i
,
t
)
∂
x
i
,where u
-
b
i, v
-
b
i, and w
-
b
i are the axial and transverse displacement components with respect to the C
i
x
i
y
i
z
i system of axes; u
b
i
=
u
b
i
(
x
i
), v
b
i
=
v
b
i
(
x
i
), and w
b
i
=
w
b
i
(
x
i
) are the corresponding components of the centroid C
i; θ
b
y
i
=
θ
b
y
i
(
x
i
) and θ
b
z
i
=
θ
b
z
i
(
x
i
) are the angles of rotation of the cross-section due to bending, with respect to its centroid; ∂
θ
b
x
i
/
d
x
i denotes the rate of change of the angle of twist θ
b
x
i
(
x
i
) regarded as the torsional curvature and φ
S
P
i is the primary warping function with respect to the cross-section’s shear center (coinciding with its centroid). Employing again the strain-displacement relations of the three-dimensional elasticity for moderate displacements [44, 45], the strain components are given as(13a)
ε
x
x
=
∂
u
-
b
i
∂
x
i
+
1
2
[
(
∂
v
-
b
i
∂
x
i
)
2
+
(
∂
w
-
b
i
∂
x
i
)
2
]
,
(13b)
γ
x
z
=
∂
w
-
b
i
∂
x
i
+
∂
u
-
b
i
∂
z
i
+
(
∂
v
-
b
i
∂
x
i
∂
v
-
b
i
∂
z
i
+
∂
w
-
b
i
∂
x
i
∂
w
-
b
i
∂
z
i
)
,
(13c)
γ
x
y
=
∂
v
-
b
i
∂
x
i
+
∂
u
-
b
i
∂
y
i
+
(
∂
v
-
b
i
∂
x
i
∂
v
-
b
i
∂
y
i
+
∂
w
-
b
i
∂
x
i
∂
w
-
b
i
∂
y
i
)
,
(13d)
ε
y
y
=
ε
z
z
=
γ
y
z
=
0
.Employing the Hooke’s stress-strain law and integrating the arising stress components over the beam’s cross-section after ignoring the nonlinear terms with respect to the angle of twist and its derivatives, the stress resultants of the beam are derived as(14a)
N
b
i
=
E
b
i
A
b
i
[
∂
u
b
i
∂
x
i
+
1
2
(
(
∂
v
b
i
∂
x
i
)
2
+
(
∂
w
b
i
∂
x
i
)
2
)
]
,
(14b)
M
b
y
i
=
-
E
b
i
I
y
i
∂
2
w
b
i
∂
x
i
2
,
(14c)
M
b
z
i
=
E
b
i
I
z
i
∂
2
v
b
i
∂
x
i
2
,
(14d)
M
b
t
P
i
=
G
b
i
I
t
i
∂
θ
b
x
i
∂
x
i
,
(14e)
M
b
w
i
=
-
E
b
i
C
S
i
∂
2
θ
b
x
i
∂
x
i
2
,where M
b
t
P
i is the primary twisting moment [15] resulting from the primary shear stress distribution; M
b
w
i is the warping moment due to torsional curvature. Furthermore I
y
i and I
z
i are the principal moments of inertia; I
S
i is the polar moment of inertia, while I
t
i and C
S
i are the torsion and warping constants of the ith beam with respect to the cross-section’s shear center (coinciding with its centroid), respectively, given as [46] (15a)
I
t
i
=
∫
Ω
(
y
i
2
+
z
i
2
+
y
i
∂
φ
S
P
i
∂
z
i
-
z
i
∂
φ
S
P
i
∂
y
i
)
d
Ω
,
(15b)
C
S
i
=
∫
Ω
(
φ
S
P
i
)
2
d
Ω
.On the basis of Hamilton’s principle, the differential equations of motion in terms of displacements are obtained as(16a)
-
∂
N
b
i
∂
x
i
+
ρ
b
i
Α
b
i
u
¨
b
i
=
∑
j
=
1
2
q
x
j
i
,
(16b)
-
N
b
i
∂
2
v
b
i
∂
x
i
2
+
∂
2
M
b
z
i
∂
x
i
2
+
ρ
b
i
Α
b
i
v
¨
b
i
-
ρ
b
i
I
z
i
∂
2
v
¨
b
i
∂
x
2
-
ρ
b
i
Α
b
i
u
¨
b
i
∂
v
b
i
∂
x
i
=
∑
j
=
1
2
(
q
y
j
i
-
q
x
j
i
∂
v
b
i
∂
x
i
-
∂
m
b
z
j
i
∂
x
i
)
,
(16c)
-
N
b
i
∂
2
w
b
i
∂
x
i
2
-
∂
2
M
b
y
i
∂
x
i
2
+
ρ
b
i
Α
b
i
w
¨
b
i
-
ρ
b
i
I
y
i
∂
2
w
¨
b
i
∂
x
2
-
ρ
b
i
Α
b
i
u
¨
b
i
∂
w
b
i
∂
x
i
=
∑
j
=
1
2
(
q
z
j
i
-
q
x
j
i
∂
w
b
i
∂
x
i
+
∂
m
b
y
j
i
∂
x
i
)
,
(16d)
-
∂
M
b
t
i
∂
x
i
-
∂
2
M
b
w
i
∂
x
i
2
+
ρ
b
i
I
S
i
θ
¨
b
x
i
-
ρ
b
i
C
S
i
∂
2
θ
¨
b
x
i
∂
x
2
=
∑
j
=
1
2
[
m
b
x
j
i
+
∂
m
b
w
j
i
∂
x
i
]
.Substituting the expressions of the stress resultants of (14a)–(14e) in (16a)–(16d), the differential equations of motion are obtained as(17a)
-
E
b
i
A
b
i
(
∂
2
u
b
i
∂
x
i
2
+
∂
w
b
i
∂
x
i
∂
2
w
b
i
∂
x
i
2
+
∂
v
b
i
∂
x
i
∂
2
v
b
i
∂
x
i
2
)
+
ρ
b
i
Α
b
i
u
¨
b
i
=
∑
j
=
1
2
q
x
j
i
,
(17b)
E
b
i
I
z
i
∂
4
v
b
i
∂
x
i
4
-
N
b
i
∂
2
v
b
i
∂
x
i
2
-
ρ
b
i
I
z
i
∂
2
v
¨
b
∂
x
2
+
ρ
b
i
Α
b
i
v
¨
b
-
ρ
b
i
Α
b
i
u
¨
b
i
∂
v
b
i
∂
x
i
=
∑
j
=
1
2
(
q
y
j
i
-
q
x
j
i
∂
v
b
i
∂
x
i
-
∂
m
b
z
j
i
∂
x
i
)
,
(17c)
E
b
i
I
y
i
∂
4
w
b
i
∂
x
i
4
-
N
b
i
∂
2
w
b
i
∂
x
i
2
-
ρ
b
i
I
y
i
∂
2
w
¨
b
∂
x
2
+
ρ
b
i
Α
b
i
w
¨
b
-
ρ
b
i
Α
b
i
u
¨
b
i
∂
w
b
i
∂
x
i
=
∑
j
=
1
2
(
q
z
j
i
-
q
x
j
i
∂
w
b
i
∂
x
i
+
∂
m
b
y
j
i
∂
x
i
)
,
(17d)
E
b
i
C
S
i
∂
4
θ
b
x
i
∂
x
i
4
-
G
b
i
I
t
i
∂
2
θ
b
x
i
∂
x
i
2
+
ρ
b
i
I
S
i
θ
¨
b
x
i
-
ρ
b
i
C
S
i
∂
2
θ
¨
b
x
i
∂
x
2
=
∑
j
=
1
2
[
m
b
x
j
i
+
∂
m
b
w
j
i
∂
x
i
]
.Moreover, the corresponding boundary conditions of the ith beam at its ends x
i
=
0
,
l
i are given as
(18)
a
b
1
i
u
b
i
+
α
b
2
i
N
b
i
=
α
b
3
i
,
(19a)
β
b
1
i
v
b
i
+
β
b
2
i
R
b
y
i
=
β
b
3
i
,
(19b)
β
-
b
1
i
θ
b
z
i
+
β
-
b
2
i
M
b
z
i
=
β
-
b
3
i
,
(20a)
γ
b
1
i
w
b
i
+
γ
b
2
i
R
b
z
i
=
γ
b
3
i
,
(20b)
γ
-
b
1
i
θ
b
y
i
+
γ
-
b
2
i
M
b
y
i
=
γ
-
b
3
i
,
(21a)
δ
b
1
i
θ
b
x
i
+
δ
b
2
i
M
b
t
i
=
δ
b
3
i
,
(21b)
δ
-
b
1
i
d
θ
b
x
i
d
x
i
+
δ
-
b
2
i
M
b
w
i
=
δ
-
b
3
iand the initial conditions as(22a)
w
b
i
(
x
i
,
0
)
=
w
b
0
i
(
x
i
)
,
(22b)
w
˙
b
i
(
x
i
,
0
)
=
w
-
b
0
i
(
x
i
)
,where the angles of rotation of the cross-section due to bending θ
b
y
i, θ
b
z
i are given from (12d) and (12e), R
b
y
i, R
b
z
i and M
b
z
i, M
b
y
i are the reactions and bending moments with respect to y
i, z
i axes, respectively, which after applying the aforementioned simplifications are given as(23a)
R
b
y
i
=
N
b
i
∂
v
b
i
∂
x
i
-
E
b
i
I
z
i
∂
3
v
b
i
∂
x
i
3
,
(23b)
R
b
z
i
=
N
b
i
∂
w
b
i
∂
x
i
-
E
b
i
I
y
i
∂
3
w
b
i
∂
x
i
3
,
(24a)
M
b
z
i
=
E
b
i
I
z
i
∂
2
v
b
i
∂
x
i
2
,
(24b)
M
b
y
i
=
-
E
b
i
I
y
i
∂
2
w
b
i
∂
x
i
2
,and M
b
t
i, M
b
w
i are the torsional and warping moments at the boundaries of the beam, respectively, given as(25a)
M
b
t
i
=
G
b
i
I
t
i
∂
θ
b
x
i
∂
x
i
-
E
b
i
C
S
i
∂
3
θ
b
x
i
∂
x
i
3
,
(25b)
M
b
w
i
=
-
E
b
i
C
S
i
∂
2
θ
b
x
i
∂
x
i
2
.Finally, α
b
k
i
,
β
b
k
i
,
β
-
b
k
i
,
γ
b
k
i
,
γ
-
b
k
i
,
δ
b
k
i
,
δ
-
b
k
i (k
=
1,2
,
3) are functions specified at the ith beam ends (x
i
=
0
,
l
i). The boundary conditions (18)–(21b) are the most general boundary conditions for the beam problem including also the elastic support. It is apparent that all types of the conventional boundary conditions (clamped, simply supported, free, or guided edge) can be derived from these equations by specifying appropriately the aforementioned coefficients.Equations (7a)–(17d) constitute a set of seven coupled and nonlinear partial differential equations including thirteen unknowns, namely u
p, v
p, w
p, u
b
i, v
b
i, w
b
i, θ
b
x
i, q
x
1
i, q
y
1
i, q
z
1
i
,
q
x
2
i, q
y
2
i, and q
z
2
i. Six additional equations are required, which result from the displacement continuity conditions in the directions of x
i, y
i, and z
i local axes along the two interface lines of each (ith) plate-beam interface. Taking into account the displacement fields expressed by (1a)–(1c) and (12a)–(12e) the displacement continuity conditions [47] can be expressed as follows.In the direction ofx
i local axis: (26a)
u
p
1
i
-
u
b
i
=
h
p
2
∂
w
p
1
i
∂
x
+
h
b
i
2
∂
w
b
i
∂
x
i
+
b
f
i
4
∂
v
b
i
∂
x
i
+
∂
θ
b
x
i
∂
x
i
(
φ
S
P
i
)
f
1
+
q
x
1
i
k
x
1
i
[
1
-
1
2
(
∂
w
b
i
∂
x
i
)
2
]
along
interface
line
1
(
f
j
=
1
i
)
,
(26b)
u
p
2
i
-
u
b
i
=
h
p
2
∂
w
p
2
i
∂
x
+
h
b
i
2
∂
w
b
i
∂
x
i
-
b
f
i
4
∂
v
b
i
∂
x
i
+
∂
θ
b
x
i
∂
x
i
(
φ
S
P
i
)
f
2
+
q
x
2
i
k
x
2
i
[
1
-
1
2
(
∂
w
b
i
∂
x
i
)
2
]
along
interface
line
2
(
f
j
=
2
i
)
.In the direction of y
i local axis: (27a)
v
p
1
i
-
v
b
i
=
h
p
2
∂
w
p
1
i
∂
y
+
h
b
i
2
θ
b
x
i
+
q
y
1
i
k
y
1
i
[
1
-
1
2
(
θ
b
x
i
)
2
]
along
interface
line
1
(
f
j
=
1
i
)
,
(27b)
v
p
2
i
-
v
b
i
=
h
p
2
∂
w
p
2
i
∂
y
+
h
b
i
2
θ
b
x
i
+
q
y
2
i
k
y
2
i
[
1
-
1
2
(
θ
b
x
i
)
2
]
along
interface
line
2
(
f
j
=
2
i
)
.In the direction of z
i local axis: (28a)
w
p
1
i
-
w
b
i
=
-
b
f
i
4
θ
b
x
i
along
interface
line
1
(
f
j
=
1
i
)
,
(28b)
w
p
2
i
-
w
b
i
=
b
f
i
4
θ
b
x
i
along
interface
line
2
(
f
j
=
2
i
)
,where (
φ
S
P
i
)
f
j is the value of the primary warping function with respect to the shear center of the beam cross-section (coinciding with its centroid) at the point of the jth interface line of the ith plate-beam interface f
j
i and k
x
j
i, k
y
j
i are the stiffness of the arbitrarily distributed shear connectors along the x
i and y
i directions, respectively. It is noted that k
x
j
i
=
k
x
j
i
(
s
x
j
i
) and k
y
j
i
=
k
y
j
i
(
s
y
j
i
) can represent any linear or nonlinear relationship between the inplane interface forces and the interface slip s
j
i in the corresponding direction. In all of the aforementioned equations the values of the primary warping function φ
S
P
i
(
y
i
,
z
i
) should be set having the appropriate algebraic sign corresponding to the local beam axes.


## 3. Integral Representations-Numerical Solution

The solution of the presented dynamic problem requires the integration of the set of (7a)–(7c) and (17a)–(17d) subjected to the prescribed boundary and initial conditions. Moreover, the displacement continuity conditions should also be fulfilled. Due to the nonlinear and coupling character of the equations of motion, an analytical solution is out of question. Therefore, a numerical solution is derived employing the analog equation method [41], a BEM-based method. Contrary to previous research efforts where the numerical analysis is based on BEM using a lumped mass assumption model after evaluating the flexibility matrix at the mass nodal points [15], in this work, a distributed mass model is employed.In the following sections, the plate and beam problems represented by (7a)–(7c) and (17a)–(17d), respectively, are examined independently and the connection between these problems is achieved employing the displacement continuity conditions.


### 3.1. For the Plate Displacement Componentsu
p, v
p, w
p

According to the precedent analysis, the large deflection analysis of the plate becomes equivalent to establishing the inplane displacement componentsu
p and v
p having continuous partial derivatives up to the second order with respect to x, y, and the deflection w
p having continuous partial derivatives up to the fourth order with respect to x, y, and all displacement components having derivatives up to the second order with respect to t. Moreover, these displacement components must satisfy the boundary value problem described by the nonlinear and coupled governing differential equations of equilibrium (equations (7a)–(7c)) inside the domain, the conditions (equations (8a)–(8e)) at the boundary Γ and the initial conditions (equations (9a)-(9b)). Equations (7a)–(7c) and (8a)–(8e) are solved using the analog equation method [41]. More specifically, setting as u
p
1
=
u
p, u
p
2
=
v
p, u
p
3
=
w
p and applying the Laplacian operator to u
p
1, u
p
2 and the biharmonic operator to u
p
3 yields(29a)
∇
2
u
p
i
=
p
p
i
(
x
,
y
,
t
)
,
(
i
=
1,2
)
,
(29b)
∇
4
u
p
3
=
p
p
3
(
x
,
y
,
t
)
.Equations (29a) and (29b) are called analog equations and they indicate that the solution of (7a)–(7c) and (8a)–(8e) can be established by solving (29a) and (29b) under the same boundary conditions (equations (8a)–(8e)), provided that the fictitious load distributions p
p
i
(
x
,
y
)
(
i
=
1
,
2,3
) are first established. These distributions can be determined using BEM. Following the procedure presented in [47], the unknown boundary quantities u
p
i
(
P
,
t
)
,
u
p
i
,
x
(
P
,
t
)
,
u
p
i
,
x
x
(
P
,
t
), and u
p
i
,
x
x
x
(
P
,
t
) (P
∈
boundary, i
=
1
,
2,3) can be expressed in terms of p
p
i (i
=
1,2
,
3) after applying the integral representations of the displacement components u
i
(
i
=
1
,
2,3
) and their derivatives with respect to x
,
y to the boundary of the plate.By discretizing the boundary of the plate intoN boundary elements and the domain Ω into M domain cells and employing the constant element assumption (as the numerical implementation becomes very simple and the obtained results are of high accuracy), the application of the integral representations and the corresponding one of the Laplacian ∇
2
u
p
3 and of the boundary conditions (8a)–(8e) to the N boundary nodal points results in a set of 8
×
N nonlinear algebraic equations relating the unknown boundary quantities with the fictitious load distributions p
p
i (i
=
1,2
,
3) that can be written as
(30)
[
E
p
11
0
0
0
E
p
22
0
0
0
E
p
33
]
{
d
p
1
d
p
2
d
p
3
}
+
{
0
D
p
1
n
l
0
D
p
2
n
l
0
0
D
p
3
n
l
}
=
{
0
α
p
3
0
β
p
3
0
0
γ
p
3
δ
p
3
}
,
where(31a)
E
p
11
=
[
A
^
1
H
1
H
^
2
0
D
p
22
D
p
23
]
,
(31b)
E
p
22
=
[
A
^
1
H
1
H
^
2
0
D
p
44
D
p
45
]
,
(31c)
E
p
33
=
[
A
2
H
1
H
2
G
1
G
2
A
1
0
0
H
1
H
2
0
D
p
78
D
p
79
0
D
p
711
0
D
p
88
D
p
89
D
p
810
0
]
.
D
p
22 to D
p
810 are N
×
N rectangular known matrices including the values of the functions a
p
j
,
β
p
j
,
γ
p
j
,
δ
p
j (j
=
1,2) of (8a)–(8d); α
p
3, β
p
3, γ
p
3, and δ
p
3 are N
×
1 known column matrices including the boundary values of the functions a
p
3
,
β
p
3
,
γ
p
3
,
δ
p
3, of (8a)–(8d); H
i (i
=
1,2), H
^
2, and G
i (i
=
1,2) are rectangular N
×
N known coefficient matrices resulting from the values of kernels at the boundary elements of the plate; A
1, A
^
1, and A
2 are N
×
M rectangular known matrices originating from the integration of kernels on the domain cells of the plate; D
p
i
n
l (i
=
1,2) are N
×
1 and D
p
3
n
l is 2
N
×
1 column matrices containing the nonlinear terms included in the expressions of the boundary conditions (equations (8a)–(8d)). It is noted that the derivatives of the unknown boundary quantities with respect to the arc length s, appearing in (8a)–(8d), are approximated employing appropriate central, backward, or forward finite difference schemes. Finally(32a)
d
p
i
=
{
p
p
i
u
^
p
i
u
^
p
i
,
n
}
T
,
(
i
=
1,2
)
,
(32b)
d
p
3
=
{
p
p
3
u
^
p
3
u
^
p
3
,
n
∇
2
u
^
p
3
(
∇
2
u
^
p
3
)
,
n
}
Tare generalized unknown vectors, where (33a)
u
^
p
i
=
{
(
u
p
i
)
1
(
u
p
i
)
2
⋯
(
u
p
i
)
N
}
T
,
(
i
=
1
,
2,3
)
(33b)
u
^
p
i
,
n
=
{
(
∂
u
p
i
∂
n
)
1
(
∂
u
p
i
∂
n
)
2
⋯
(
∂
u
p
i
∂
n
)
N
}
T
,
(
i
=
1
,
2,3
)
(33c)
∇
2
u
^
p
3
=
{
(
∇
2
u
p
3
)
1
(
∇
2
u
p
3
)
2
⋯
(
∇
2
u
p
3
)
N
}
T
,
(33d)
(
∇
2
u
^
p
3
)
,
n
=
{
(
∂
∇
2
u
p
3
∂
n
)
1
(
∂
∇
2
u
p
3
∂
n
)
2
⋯
(
∂
∇
2
u
p
3
∂
n
)
N
}
Tare vectors including the unknown boundary values of the respective boundary quantities and p
p
i
=
{
(
p
p
i
)
1
(
p
p
i
)
2
⋯
(
p
p
i
)
M
}
T, (i
=
1
,
2,3), are vectors containing the M unknown nodal values of the fictitious loads at the domain cells of the plate. In the case that the boundary Γ has k free or transversely elastically restrained corners, k additional equations must be satisfied together with (30). These additional equations result from the application of the corner condition (8e) on the k corners following the procedure presented in [48].Discretization of the integral representations of the displacement componentsu
i and their derivatives with respect to x, y [49] and application to the M domain nodal points yields(34a)
u
p
i
=
B
p
i
d
p
i
,
(
i
=
1,2
,
3
)
,
(34b)
u
p
i
,
x
=
B
p
i
,
x
d
p
i
,
(
i
=
1,2
,
3
)
,
(34c)
u
p
i
,
y
=
B
p
i
,
y
d
p
i
,
(
i
=
1,2
,
3
)
,
(34d)
u
p
i
,
x
x
=
B
p
i
,
x
x
d
p
i
,
(
i
=
1,2
,
3
)
,
(34e)
u
p
i
,
y
y
=
B
p
i
,
y
y
d
p
i
,
(
i
=
1,2
,
3
)
,
(34f)
u
p
i
,
x
y
=
B
p
i
,
x
y
d
p
i
,
(
i
=
1,2
,
3
)
,where B
p
i, B
p
i
,
x, B
p
i
,
y, B
p
i
,
x
x, B
p
i
,
y
y, B
p
i
,
x
y (i
=
1,2) are M
×
(
2
N
+
M
) and B
p
3, B
p
3
,
x, B
p
3
,
y, B
p
3
,
x
x, B
p
3
,
y
y, B
p
3
,
x
y are M
×
(
4
N
+
M
) known coefficient matrices. In evaluating the domain integrals of the kernels over the domain cells, singular and hypersingular integrals arise. They are computed by transforming them to line integrals on the boundary of the cell [50, 51]. The final step of the AEM is to apply (7a)–(7c) to the M nodal points inside Ω, yielding the following equations:(35a)
G
p
h
p
[
p
p
1
+
1
+
ν
p
1
-
ν
p
(
B
p
1
,
x
x
d
p
1
+
B
p
2
,
x
y
d
p
2
)
h
h
h
h
h
h
+
{
2
1
-
ν
p
(
B
p
3
,
x
x
d
p
3
)
d
g
B
p
3
,
x
d
p
3
h
h
h
h
h
h
h
h
+
(
B
p
3
,
y
y
d
p
3
)
d
g
·
B
p
3
,
x
d
p
3
2
1
-
ν
p
}
h
h
h
h
h
h
+
1
+
ν
p
1
-
ν
p
(
B
p
3
,
x
y
d
p
3
)
d
g
B
p
3
,
y
d
p
3
]
-
ρ
p
h
p
B
p
1
d
¨
p
1
=
Z
q
x
,
(35b)
G
p
h
p
[
p
p
2
+
1
+
ν
p
1
-
ν
p
(
B
p
1
,
x
y
d
p
1
+
B
p
2
,
y
y
d
p
2
)
h
h
h
h
h
h
+
{
2
1
-
ν
p
(
B
p
3
,
y
y
d
p
3
)
d
g
B
p
3
,
y
d
p
3
h
h
h
h
h
h
h
h
h
+
(
B
p
3
,
x
x
d
p
3
)
d
g
·
B
p
3
,
y
d
p
3
2
1
-
ν
p
}
h
h
h
h
h
h
+
1
+
ν
p
1
-
ν
p
(
B
p
3
,
x
y
d
p
3
)
d
g
B
p
3
,
x
d
p
3
]
-
ρ
p
h
p
B
p
2
d
¨
p
2
=
Z
q
y
,
(35c)
D
p
p
3
-
C
×
{
1
2
[
(
B
p
1
,
x
d
p
1
)
d
g
B
p
3
,
x
x
d
p
3
hhhhhhhh
+
1
2
(
B
p
3
,
x
d
p
3
)
d
g
(
B
p
3
,
x
d
p
3
)
d
g
B
p
3
,
x
x
d
p
3
hhhhhhhh
+
ν
p
(
B
p
2
,
y
d
p
2
)
d
g
B
p
3
,
x
x
d
p
3
+
1
2
ν
p
(
B
p
3
,
y
d
p
3
)
d
g
hhhhhhhh
×
(
B
p
3
,
y
d
p
3
)
d
g
B
p
3
,
x
x
d
p
3
]
hhhhhhh
+
(
1
-
ν
p
)
hhhhhhh
×
[
(
B
p
1
,
y
d
p
1
)
d
g
B
p
3
,
x
y
d
p
3
hhhhhhhhhhh
+
(
B
p
2
,
x
d
p
2
)
d
g
B
p
3
,
x
y
d
p
3
hhhhhhhhhhh
+
(
B
p
3
,
x
d
p
3
)
d
g
·
(
B
p
3
,
y
d
p
3
)
d
g
B
p
3
,
x
y
d
p
3
]
hhhhhhhh
+
[
(
B
p
2
,
y
d
p
2
)
d
g
B
p
3
,
y
y
d
p
3
+
1
2
(
B
p
3
,
y
d
p
3
)
d
g
hhhhhhhhhhh
·
(
B
p
3
,
y
d
p
3
)
d
g
B
p
3
,
y
y
d
p
3
hhhhhhhhhhh
+
ν
p
(
B
p
1
,
x
d
p
1
)
d
g
B
p
3
,
y
y
d
p
3
hhhhhhhhhhh
+
1
2
ν
p
(
B
p
3
,
x
d
p
3
)
d
g
hhhhhhhhhhh
·
(
B
p
3
,
x
d
p
3
)
d
g
B
p
3
,
y
y
d
p
3
1
2
]
}
+
ρ
p
h
p
B
p
3
d
¨
p
3
-
ρ
p
h
p
(
B
p
1
d
¨
p
1
)
d
g
B
p
3
,
x
d
p
3
-
ρ
p
h
p
(
B
p
2
d
¨
p
2
)
d
g
B
p
3
,
y
d
p
3
-
ρ
p
h
p
3
12
B
p
3
,
x
x
d
p
3
-
ρ
p
h
p
3
12
B
p
3
,
y
y
d
p
3
=
g
-
Z
q
z
-
Z
X
y
q
y
-
Z
X
x
q
x
+
(
Z
q
x
)
d
g
B
p
3
,
x
d
p
3
+
(
Z
q
y
)
d
g
B
p
3
,
y
d
p
3
,where q
x
=
{
q
x
1
q
x
2
}
T, {
q
y
}
=
{
q
y
1
q
y
2
}
T, and {
q
z
}
=
{
q
z
1
q
z
2
}
T are vectors of dimension 2L including the unknown q
x
j
i, q
y
j
i, q
z
j
i (j
=
1,2) interface forces; 2
L is the total number of the nodal points at each plate-beam interface; Z is a position M
×
2
L matrix which converts the vectors q
x, q
y, q
z into corresponding ones with length M; the symbol (
·
)
d
g indicates a diagonal M
×
M matrix with the elements of the included column matrix. Matrices X
x and X
y of dimension 2
L
×
2
L result after approximating the bending moment derivatives of m
p
y
j
i
=
q
x
j
i
h
p
/
2, m
p
x
j
i
=
q
y
j
i
h
p
/
2, respectively, using appropriately central, backward, or forward differences.


### 3.2. For the Beam Displacement Componentsu
b
i, v
b
i, and w
b
i and for the Angle of Twist θ
b
x
i

According to the precedent analysis, the nonlinear dynamic analysis of theith beam reduces in establishing the displacement components u
b
i
(
x
i
,
t
) and v
b
i
(
x
i
,
t
), w
b
i
(
x
i
,
t
), θ
b
x
i
(
x
i
,
t
) having continuous derivatives up to the second order and up to the fourth order with respect to x, respectively, and also having derivatives up to the second order with respect to t. Moreover, these displacement components must satisfy the coupled governing differential equations (17a)–(17d) inside the beam, the boundary conditions at the beam ends (18)–(21b), and the initial conditions (22a) and (22b).Letu
b
1
i
=
u
b
i
(
x
i
,
t
), u
b
2
i
=
v
b
i
(
x
i
,
t
), u
b
3
i
=
w
b
i
(
x
i
,
t
), and u
b
4
i
=
θ
b
x
i
(
x
i
,
t
) be the sought solutions of the problem represented by (17a)–(17d). Differentiating with respect to x these functions, two and four times, respectively, yields(36a)
∂
2
u
b
1
i
∂
x
i
2
=
p
b
1
i
(
x
i
,
t
)
,
(36b)
∂
4
u
b
j
i
∂
x
i
4
=
p
b
j
i
(
x
i
,
t
)
,
(
j
=
2,3
,
4
)
.Equations (36a) and (36b) are quasi-static and indicate that the solution of (17a)–(17d) can be established by solving (36a) and (36b) under the same boundary conditions (equations (18)–(21b)) provided that the fictitious load distributions p
b
j
i
(
x
i
,
t
)
(
j
=
1
,
2,3
,
4
) are first established. These distributions can be determined following the procedure presented in [49] and employing the constant element assumption for the load distributions p
b
j
i along the L internal beam elements (as the numerical implementation becomes very simple and the obtained results are of high accuracy). Thus, the integral representations of the displacement components u
b
j
i (j
=
1
,
2,3
,
4) and their derivatives with respect to x
i when applied to the beam ends (0
,
l
i), together with the boundary conditions (18)–(21b), are employed to express the unknown boundary quantities u
b
j
i
(
ζ
i
), u
b
j
,
x
i
(
ζ
i
), u
b
j
,
x
x
i
(
ζ
i
), and u
b
j
,
x
x
x
i
(
ζ
i
) (ζ
i
=
0
,
l
i) in terms of p
b
j
i (j
=
1
,
2,3
,
4). Thus, the following set of 28 nonlinear algebraic equations for the ith beam is obtained as
(37)
[
E
b
11
i
0
0
0
0
E
b
22
i
0
0
0
0
E
b
33
i
0
0
0
0
E
b
44
i
]
{
d
b
1
i
d
b
2
i
d
b
3
i
d
b
4
i
}
+
{
0
D
b
1
i
n
l
0
0
D
b
2
i
n
l
0
0
D
b
3
i
n
l
0
0
0
}
=
{
0
α
b
3
i
0
0
β
b
3
i
0
0
γ
b
3
i
0
0
δ
b
3
i
}
,
where E
b
11
i and E
b
22
i
,
E
b
33
i
,
E
b
44
i are known matrices of dimension 4
×
(
N
+
4
) and 8
×
(
N
+
8
), respectively; D
b
1
i
n
l is a 2
×
1 and D
b
j
i
n
l (j
=
2,3) are 4
×
1 column matrices containing the nonlinear terms included in the expressions of the boundary conditions (equations (18)–(20b)); α
b
3
i and β
b
3
i, γ
b
3
i, δ
b
3
i are 2
×
1 and 4
×
1 known column matrices, respectively, including the boundary values of the functions a
b
3
i and β
b
3
i
,
β
-
b
3
i
,
γ
b
3
i
,
γ
-
b
3
i
,
δ
b
3
i
,
δ
-
b
3
i of (18)–(21b). Finally(38a)
d
b
1
i
=
{
p
b
1
i
u
^
b
1
i
u
^
b
1
,
x
i
}
T
,
(38b)
d
b
j
i
=
{
p
b
j
i
u
^
b
j
i
u
^
b
j
,
x
i
u
^
b
j
,
x
x
i
u
^
b
j
,
x
x
x
i
}
T
,
(
j
=
2,3
,
4
)are generalized unknown vectors, where(39a)
u
^
b
j
i
=
{
u
b
j
i
(
0
,
t
)
u
b
j
i
(
l
i
,
t
)
}
T
,
(
j
=
1
,
2,3
,
4
)
,
(39b)
u
^
b
j
,
x
i
=
{
∂
u
b
j
i
(
0
,
t
)
∂
x
i
∂
u
b
j
i
(
l
i
,
t
)
∂
x
i
}
T
,
(
j
=
1
,
2,3
,
4
)
,
(39c)
u
^
b
j
,
x
x
i
=
{
∂
2
u
b
j
i
(
0
,
t
)
∂
x
i
2
∂
2
u
b
j
i
(
l
i
,
t
)
∂
x
i
2
}
T
,
(
j
=
2,3
,
4
)
,
(39d)
u
^
b
j
,
x
x
x
i
=
{
∂
3
u
b
j
i
(
0
,
t
)
∂
x
i
3
∂
3
u
b
j
i
(
l
i
,
t
)
∂
x
i
3
}
T
,
(
j
=
2,3
,
4
)are vectors including the two unknown boundary values of the respective boundary quantities and p
b
j
i
=
{
(
p
b
j
i
)
1
(
p
b
j
i
)
2
⋯
(
p
b
j
i
)
L
}
T
(
j
=
1
,
2,3
,
4
) are vectors including the L unknown nodal values of the fictitious loads.Discretization of the integral representations of the unknown quantitiesu
b
j
i (j
=
1
,
2,3
,
4) and those of their derivatives with respect to x
i inside the beam x
i
∈
(
0
,
l
i
) and application to the L collocation nodal points yields(40a)
u
b
1
i
=
B
b
1
i
d
b
1
i
,
(40b)
u
b
1
,
x
i
=
B
b
1
,
x
i
d
b
1
i
,
(40c)
u
b
j
i
=
B
b
j
i
d
b
j
i
,
(
j
=
2,3
,
4
)
,
(40d)
u
b
j
,
x
i
=
B
b
j
,
x
i
d
b
j
i
,
(
j
=
2,3
,
4
)
,
(40e)
u
b
j
,
x
x
i
=
B
b
j
,
x
x
i
d
b
j
i
,
(
j
=
2,3
,
4
)
,
(40f)
u
b
j
,
x
x
x
i
=
B
b
j
,
x
x
x
i
d
b
j
i
,
(
j
=
2,3
,
4
)
,where B
b
1
i
,
B
b
1
,
x
i are L
×
(
L
+
4
) and B
b
j
i
,
B
b
j
,
x
i
,
B
b
j
,
x
x
i
,
B
b
j
,
x
x
x
i (j
=
2,3
,
4) are L
×
(
L
+
8
) known coefficient matrices, respectively.Applying (17a)–(17d) to the L collocation points and employing (40a)–(40f), 4
×
L nonlinear algebraic equations for each (ith) beam are formulated as(41a)
-
E
b
i
A
b
i
[
p
b
1
i
+
(
B
b
2
,
x
i
d
b
2
i
)
d
g
B
b
2
,
x
x
i
d
b
2
i
hhhhhhh
+
(
B
b
3
,
x
i
d
b
3
i
)
d
g
B
b
3
,
x
x
i
d
b
3
i
]
+
ρ
b
i
Α
b
i
T
1
q
¨
1
=
q
x
1
i
+
q
x
2
i
,
(41b)
E
b
i
I
z
i
p
b
2
i
-
(
N
b
i
)
d
g
B
b
2
,
x
x
i
d
b
2
i
-
ρ
b
i
I
z
i
B
b
2
,
x
x
i
d
¨
b
2
i
+
ρ
b
i
Α
b
i
B
b
2
i
d
¨
b
2
i
-
ρ
b
i
Α
b
i
(
B
b
1
i
d
¨
b
1
i
)
d
g
B
b
2
,
x
i
d
b
2
i
=
q
y
1
i
+
q
y
2
i
-
(
B
b
2
,
x
i
d
b
2
i
)
d
g
(
q
x
1
i
+
q
x
2
i
)
-
X
b
y
i
(
q
x
1
i
+
q
x
2
i
)
,
(41c)
E
b
i
I
y
i
p
b
3
i
-
(
N
b
i
)
d
g
B
b
3
,
x
x
i
d
b
3
i
-
ρ
b
i
I
y
i
B
b
3
,
x
x
i
d
¨
b
3
i
+
ρ
b
i
Α
b
i
B
b
3
i
d
¨
b
3
i
-
ρ
b
i
Α
b
i
(
B
b
1
i
d
¨
b
1
i
)
d
g
B
b
3
,
x
i
d
b
3
i
=
q
z
1
i
+
q
z
2
i
-
(
B
b
3
,
x
i
d
b
3
i
)
d
g
(
q
x
1
i
+
q
x
2
i
)
+
X
b
z
i
(
q
x
1
i
+
q
x
2
i
)
,
(41d)
E
b
i
C
S
i
p
b
4
i
-
G
b
i
I
t
i
B
b
4
,
x
x
i
d
b
4
i
+
ρ
b
i
I
p
i
B
b
4
i
d
¨
b
4
i
-
ρ
b
i
C
S
i
B
b
4
,
x
x
i
d
¨
b
4
i
=
e
y
1
i
q
z
1
i
+
e
y
2
i
q
z
2
i
-
e
z
1
i
q
y
1
i
-
e
z
2
i
q
y
2
i
+
Χ
b
w
i
(
q
x
1
i
+
q
x
2
i
)
,where (
N
b
i
)
d
g is a diagonal L
×
L matrix including the values of the axial forces of the ith beam; the symbol (
·
)
d
g indicates a diagonal L
×
L matrix with the elements of the included column matrix. The matrices X
b
z
i, X
b
y
i, X
b
w
i result after approximating the derivatives of m
b
y
j
i, m
b
z
j
i, m
b
w
j
i using appropriately central, backward, or forward differences. Their dimensions are also L
×
L. Moreover e
y
1
i, e
y
2
i, e
z
1
i, e
z
2
i are diagonal L
×
L matrices including the values of the eccentricities e
y
j
i, e
z
j
i of the components q
z
j
i, q
y
j
i with respect to the ith beam shear center axis (coinciding with its centroid), respectively.Employing (34a)–(34f) and (40a)–(40f), the discretized counterpart of (26a)-(26b), (27a)-(27b) and (28a)-(28b) at the L nodal points of each interface, is written as(42a)
Y
1
B
p
1
d
p
1
-
B
b
1
i
d
b
1
i
=
h
p
2
Y
1
B
p
3
,
x
d
p
3
+
h
b
i
2
B
b
3
,
x
i
d
b
3
+
b
f
i
4
B
b
2
,
x
i
d
b
2
+
(
φ
S
P
i
)
f
1
B
b
4
,
x
i
d
b
4
+
K
x
1
i
[
1
-
1
2
(
B
b
3
,
x
i
d
b
3
)
d
g
(
B
b
3
,
x
i
d
b
3
)
d
g
]
(
q
x
1
i
+
q
x
2
i
)
,
(42b)
Y
2
B
p
1
d
p
1
-
B
b
1
i
d
b
1
i
=
h
p
2
Y
2
B
p
3
,
x
d
p
3
+
h
b
i
2
B
b
3
,
x
i
d
b
3
-
b
f
i
4
B
b
2
,
x
i
d
b
2
+
(
φ
S
P
i
)
f
2
B
b
4
,
x
i
d
b
4
+
K
x
2
i
[
1
-
1
2
(
B
b
3
,
x
i
d
b
3
)
d
g
(
B
b
3
,
x
i
d
b
3
)
d
g
]
(
q
x
1
i
+
q
x
2
i
)
,
(43a)
Y
1
B
p
2
d
p
2
-
B
b
2
i
d
b
2
i
=
h
p
2
Y
1
B
p
3
,
y
d
p
3
+
h
b
i
2
B
b
4
i
d
b
4
i
+
K
y
1
i
[
1
-
1
2
(
B
b
4
,
x
i
d
b
4
)
d
g
(
B
b
4
,
x
i
d
b
4
)
d
g
]
(
q
x
1
i
+
q
x
2
i
)
,
(43b)
Y
2
B
p
2
d
p
2
-
B
b
2
i
d
b
2
i
=
h
p
2
Y
2
B
p
3
,
y
d
p
3
+
h
b
i
2
B
b
4
i
d
b
4
i
+
K
y
2
i
[
1
-
1
2
(
B
b
4
,
x
i
d
b
4
)
d
g
(
B
b
4
,
x
i
d
b
4
)
d
g
]
(
q
x
1
i
+
q
x
2
i
)
,
(44a)
Y
1
B
p
3
d
p
3
-
B
b
3
i
d
b
3
i
=
-
b
f
i
4
B
b
4
i
d
b
4
i
,
(44b)
Y
2
B
p
3
d
p
3
-
B
b
3
i
d
b
3
i
=
b
f
i
4
B
b
4
i
d
b
4
i
.Equations (30), (35a)–(35c), (37), and (41a)–(41d) together with continuity conditions (42a), (42b), (43a), (43b), (44a), and (44b) constitute a nonlinear system of algebraic equations with respect to q
x
1, q
x
2, q
y
1, q
y
2, q
z
1, and q
z
2 (interface forces) and d
p
i (i
=
1,2
,
3) and d
b
j
i (i
=
1
⋯
I), (j
=
1,2
,
3,4) (generalized unknown vectors of the plate and the beams). This system is solved using iterative numerical methods [52]. It is worth noting here that Y
1, Y
2 are position L
×
M matrices which convert the matrices B
p
i (i
=
1,2
,
3), B
p
3
,
x, and B
p
3
,
y into corresponding ones with dimensions L
×
M, appropriately referring to the nodal points of the two interface lines f
j
=
1
i, f
j
=
2
i, respectively. Moreover, K
x
1
i, K
x
2
i, K
y
1
i, and K
y
2
i are diagonal L
×
L matrices including the value of the flexibilities 1
/
k
x
j
i and 1
/
k
y
j
i
(
j
=
1,2
) of the arbitrary distributed shear connectors along the x
i and y
i directions, respectively.Finally, it is worth noting that beams placed along the boundary of the plate are treated as every other stiffening beam, since the lines of action of the integrated interface force componentsq
x
j
i, q
y
j
i, and q
z
j
i
(
j
=
1,2
) will also be internal ones, taking special care during the numerical evaluation of the line integrals in order to avoid their “near singular integral behaviour.” According to this, boundary elements that are very close to each other (distance smaller than their length) are divided in subelements, in each of which Gauss integration is applied [51].


## 3.1. For the Plate Displacement Componentsu
p, v
p, w
p

According to the precedent analysis, the large deflection analysis of the plate becomes equivalent to establishing the inplane displacement componentsu
p and v
p having continuous partial derivatives up to the second order with respect to x, y, and the deflection w
p having continuous partial derivatives up to the fourth order with respect to x, y, and all displacement components having derivatives up to the second order with respect to t. Moreover, these displacement components must satisfy the boundary value problem described by the nonlinear and coupled governing differential equations of equilibrium (equations (7a)–(7c)) inside the domain, the conditions (equations (8a)–(8e)) at the boundary Γ and the initial conditions (equations (9a)-(9b)). Equations (7a)–(7c) and (8a)–(8e) are solved using the analog equation method [41]. More specifically, setting as u
p
1
=
u
p, u
p
2
=
v
p, u
p
3
=
w
p and applying the Laplacian operator to u
p
1, u
p
2 and the biharmonic operator to u
p
3 yields(29a)
∇
2
u
p
i
=
p
p
i
(
x
,
y
,
t
)
,
(
i
=
1,2
)
,
(29b)
∇
4
u
p
3
=
p
p
3
(
x
,
y
,
t
)
.Equations (29a) and (29b) are called analog equations and they indicate that the solution of (7a)–(7c) and (8a)–(8e) can be established by solving (29a) and (29b) under the same boundary conditions (equations (8a)–(8e)), provided that the fictitious load distributions p
p
i
(
x
,
y
)
(
i
=
1
,
2,3
) are first established. These distributions can be determined using BEM. Following the procedure presented in [47], the unknown boundary quantities u
p
i
(
P
,
t
)
,
u
p
i
,
x
(
P
,
t
)
,
u
p
i
,
x
x
(
P
,
t
), and u
p
i
,
x
x
x
(
P
,
t
) (P
∈
boundary, i
=
1
,
2,3) can be expressed in terms of p
p
i (i
=
1,2
,
3) after applying the integral representations of the displacement components u
i
(
i
=
1
,
2,3
) and their derivatives with respect to x
,
y to the boundary of the plate.By discretizing the boundary of the plate intoN boundary elements and the domain Ω into M domain cells and employing the constant element assumption (as the numerical implementation becomes very simple and the obtained results are of high accuracy), the application of the integral representations and the corresponding one of the Laplacian ∇
2
u
p
3 and of the boundary conditions (8a)–(8e) to the N boundary nodal points results in a set of 8
×
N nonlinear algebraic equations relating the unknown boundary quantities with the fictitious load distributions p
p
i (i
=
1,2
,
3) that can be written as
(30)
[
E
p
11
0
0
0
E
p
22
0
0
0
E
p
33
]
{
d
p
1
d
p
2
d
p
3
}
+
{
0
D
p
1
n
l
0
D
p
2
n
l
0
0
D
p
3
n
l
}
=
{
0
α
p
3
0
β
p
3
0
0
γ
p
3
δ
p
3
}
,
where(31a)
E
p
11
=
[
A
^
1
H
1
H
^
2
0
D
p
22
D
p
23
]
,
(31b)
E
p
22
=
[
A
^
1
H
1
H
^
2
0
D
p
44
D
p
45
]
,
(31c)
E
p
33
=
[
A
2
H
1
H
2
G
1
G
2
A
1
0
0
H
1
H
2
0
D
p
78
D
p
79
0
D
p
711
0
D
p
88
D
p
89
D
p
810
0
]
.
D
p
22 to D
p
810 are N
×
N rectangular known matrices including the values of the functions a
p
j
,
β
p
j
,
γ
p
j
,
δ
p
j (j
=
1,2) of (8a)–(8d); α
p
3, β
p
3, γ
p
3, and δ
p
3 are N
×
1 known column matrices including the boundary values of the functions a
p
3
,
β
p
3
,
γ
p
3
,
δ
p
3, of (8a)–(8d); H
i (i
=
1,2), H
^
2, and G
i (i
=
1,2) are rectangular N
×
N known coefficient matrices resulting from the values of kernels at the boundary elements of the plate; A
1, A
^
1, and A
2 are N
×
M rectangular known matrices originating from the integration of kernels on the domain cells of the plate; D
p
i
n
l (i
=
1,2) are N
×
1 and D
p
3
n
l is 2
N
×
1 column matrices containing the nonlinear terms included in the expressions of the boundary conditions (equations (8a)–(8d)). It is noted that the derivatives of the unknown boundary quantities with respect to the arc length s, appearing in (8a)–(8d), are approximated employing appropriate central, backward, or forward finite difference schemes. Finally(32a)
d
p
i
=
{
p
p
i
u
^
p
i
u
^
p
i
,
n
}
T
,
(
i
=
1,2
)
,
(32b)
d
p
3
=
{
p
p
3
u
^
p
3
u
^
p
3
,
n
∇
2
u
^
p
3
(
∇
2
u
^
p
3
)
,
n
}
Tare generalized unknown vectors, where (33a)
u
^
p
i
=
{
(
u
p
i
)
1
(
u
p
i
)
2
⋯
(
u
p
i
)
N
}
T
,
(
i
=
1
,
2,3
)
(33b)
u
^
p
i
,
n
=
{
(
∂
u
p
i
∂
n
)
1
(
∂
u
p
i
∂
n
)
2
⋯
(
∂
u
p
i
∂
n
)
N
}
T
,
(
i
=
1
,
2,3
)
(33c)
∇
2
u
^
p
3
=
{
(
∇
2
u
p
3
)
1
(
∇
2
u
p
3
)
2
⋯
(
∇
2
u
p
3
)
N
}
T
,
(33d)
(
∇
2
u
^
p
3
)
,
n
=
{
(
∂
∇
2
u
p
3
∂
n
)
1
(
∂
∇
2
u
p
3
∂
n
)
2
⋯
(
∂
∇
2
u
p
3
∂
n
)
N
}
Tare vectors including the unknown boundary values of the respective boundary quantities and p
p
i
=
{
(
p
p
i
)
1
(
p
p
i
)
2
⋯
(
p
p
i
)
M
}
T, (i
=
1
,
2,3), are vectors containing the M unknown nodal values of the fictitious loads at the domain cells of the plate. In the case that the boundary Γ has k free or transversely elastically restrained corners, k additional equations must be satisfied together with (30). These additional equations result from the application of the corner condition (8e) on the k corners following the procedure presented in [48].Discretization of the integral representations of the displacement componentsu
i and their derivatives with respect to x, y [49] and application to the M domain nodal points yields(34a)
u
p
i
=
B
p
i
d
p
i
,
(
i
=
1,2
,
3
)
,
(34b)
u
p
i
,
x
=
B
p
i
,
x
d
p
i
,
(
i
=
1,2
,
3
)
,
(34c)
u
p
i
,
y
=
B
p
i
,
y
d
p
i
,
(
i
=
1,2
,
3
)
,
(34d)
u
p
i
,
x
x
=
B
p
i
,
x
x
d
p
i
,
(
i
=
1,2
,
3
)
,
(34e)
u
p
i
,
y
y
=
B
p
i
,
y
y
d
p
i
,
(
i
=
1,2
,
3
)
,
(34f)
u
p
i
,
x
y
=
B
p
i
,
x
y
d
p
i
,
(
i
=
1,2
,
3
)
,where B
p
i, B
p
i
,
x, B
p
i
,
y, B
p
i
,
x
x, B
p
i
,
y
y, B
p
i
,
x
y (i
=
1,2) are M
×
(
2
N
+
M
) and B
p
3, B
p
3
,
x, B
p
3
,
y, B
p
3
,
x
x, B
p
3
,
y
y, B
p
3
,
x
y are M
×
(
4
N
+
M
) known coefficient matrices. In evaluating the domain integrals of the kernels over the domain cells, singular and hypersingular integrals arise. They are computed by transforming them to line integrals on the boundary of the cell [50, 51]. The final step of the AEM is to apply (7a)–(7c) to the M nodal points inside Ω, yielding the following equations:(35a)
G
p
h
p
[
p
p
1
+
1
+
ν
p
1
-
ν
p
(
B
p
1
,
x
x
d
p
1
+
B
p
2
,
x
y
d
p
2
)
h
h
h
h
h
h
+
{
2
1
-
ν
p
(
B
p
3
,
x
x
d
p
3
)
d
g
B
p
3
,
x
d
p
3
h
h
h
h
h
h
h
h
+
(
B
p
3
,
y
y
d
p
3
)
d
g
·
B
p
3
,
x
d
p
3
2
1
-
ν
p
}
h
h
h
h
h
h
+
1
+
ν
p
1
-
ν
p
(
B
p
3
,
x
y
d
p
3
)
d
g
B
p
3
,
y
d
p
3
]
-
ρ
p
h
p
B
p
1
d
¨
p
1
=
Z
q
x
,
(35b)
G
p
h
p
[
p
p
2
+
1
+
ν
p
1
-
ν
p
(
B
p
1
,
x
y
d
p
1
+
B
p
2
,
y
y
d
p
2
)
h
h
h
h
h
h
+
{
2
1
-
ν
p
(
B
p
3
,
y
y
d
p
3
)
d
g
B
p
3
,
y
d
p
3
h
h
h
h
h
h
h
h
h
+
(
B
p
3
,
x
x
d
p
3
)
d
g
·
B
p
3
,
y
d
p
3
2
1
-
ν
p
}
h
h
h
h
h
h
+
1
+
ν
p
1
-
ν
p
(
B
p
3
,
x
y
d
p
3
)
d
g
B
p
3
,
x
d
p
3
]
-
ρ
p
h
p
B
p
2
d
¨
p
2
=
Z
q
y
,
(35c)
D
p
p
3
-
C
×
{
1
2
[
(
B
p
1
,
x
d
p
1
)
d
g
B
p
3
,
x
x
d
p
3
hhhhhhhh
+
1
2
(
B
p
3
,
x
d
p
3
)
d
g
(
B
p
3
,
x
d
p
3
)
d
g
B
p
3
,
x
x
d
p
3
hhhhhhhh
+
ν
p
(
B
p
2
,
y
d
p
2
)
d
g
B
p
3
,
x
x
d
p
3
+
1
2
ν
p
(
B
p
3
,
y
d
p
3
)
d
g
hhhhhhhh
×
(
B
p
3
,
y
d
p
3
)
d
g
B
p
3
,
x
x
d
p
3
]
hhhhhhh
+
(
1
-
ν
p
)
hhhhhhh
×
[
(
B
p
1
,
y
d
p
1
)
d
g
B
p
3
,
x
y
d
p
3
hhhhhhhhhhh
+
(
B
p
2
,
x
d
p
2
)
d
g
B
p
3
,
x
y
d
p
3
hhhhhhhhhhh
+
(
B
p
3
,
x
d
p
3
)
d
g
·
(
B
p
3
,
y
d
p
3
)
d
g
B
p
3
,
x
y
d
p
3
]
hhhhhhhh
+
[
(
B
p
2
,
y
d
p
2
)
d
g
B
p
3
,
y
y
d
p
3
+
1
2
(
B
p
3
,
y
d
p
3
)
d
g
hhhhhhhhhhh
·
(
B
p
3
,
y
d
p
3
)
d
g
B
p
3
,
y
y
d
p
3
hhhhhhhhhhh
+
ν
p
(
B
p
1
,
x
d
p
1
)
d
g
B
p
3
,
y
y
d
p
3
hhhhhhhhhhh
+
1
2
ν
p
(
B
p
3
,
x
d
p
3
)
d
g
hhhhhhhhhhh
·
(
B
p
3
,
x
d
p
3
)
d
g
B
p
3
,
y
y
d
p
3
1
2
]
}
+
ρ
p
h
p
B
p
3
d
¨
p
3
-
ρ
p
h
p
(
B
p
1
d
¨
p
1
)
d
g
B
p
3
,
x
d
p
3
-
ρ
p
h
p
(
B
p
2
d
¨
p
2
)
d
g
B
p
3
,
y
d
p
3
-
ρ
p
h
p
3
12
B
p
3
,
x
x
d
p
3
-
ρ
p
h
p
3
12
B
p
3
,
y
y
d
p
3
=
g
-
Z
q
z
-
Z
X
y
q
y
-
Z
X
x
q
x
+
(
Z
q
x
)
d
g
B
p
3
,
x
d
p
3
+
(
Z
q
y
)
d
g
B
p
3
,
y
d
p
3
,where q
x
=
{
q
x
1
q
x
2
}
T, {
q
y
}
=
{
q
y
1
q
y
2
}
T, and {
q
z
}
=
{
q
z
1
q
z
2
}
T are vectors of dimension 2L including the unknown q
x
j
i, q
y
j
i, q
z
j
i (j
=
1,2) interface forces; 2
L is the total number of the nodal points at each plate-beam interface; Z is a position M
×
2
L matrix which converts the vectors q
x, q
y, q
z into corresponding ones with length M; the symbol (
·
)
d
g indicates a diagonal M
×
M matrix with the elements of the included column matrix. Matrices X
x and X
y of dimension 2
L
×
2
L result after approximating the bending moment derivatives of m
p
y
j
i
=
q
x
j
i
h
p
/
2, m
p
x
j
i
=
q
y
j
i
h
p
/
2, respectively, using appropriately central, backward, or forward differences.


## 3.2. For the Beam Displacement Componentsu
b
i, v
b
i, and w
b
i and for the Angle of Twist θ
b
x
i

According to the precedent analysis, the nonlinear dynamic analysis of theith beam reduces in establishing the displacement components u
b
i
(
x
i
,
t
) and v
b
i
(
x
i
,
t
), w
b
i
(
x
i
,
t
), θ
b
x
i
(
x
i
,
t
) having continuous derivatives up to the second order and up to the fourth order with respect to x, respectively, and also having derivatives up to the second order with respect to t. Moreover, these displacement components must satisfy the coupled governing differential equations (17a)–(17d) inside the beam, the boundary conditions at the beam ends (18)–(21b), and the initial conditions (22a) and (22b).Letu
b
1
i
=
u
b
i
(
x
i
,
t
), u
b
2
i
=
v
b
i
(
x
i
,
t
), u
b
3
i
=
w
b
i
(
x
i
,
t
), and u
b
4
i
=
θ
b
x
i
(
x
i
,
t
) be the sought solutions of the problem represented by (17a)–(17d). Differentiating with respect to x these functions, two and four times, respectively, yields(36a)
∂
2
u
b
1
i
∂
x
i
2
=
p
b
1
i
(
x
i
,
t
)
,
(36b)
∂
4
u
b
j
i
∂
x
i
4
=
p
b
j
i
(
x
i
,
t
)
,
(
j
=
2,3
,
4
)
.Equations (36a) and (36b) are quasi-static and indicate that the solution of (17a)–(17d) can be established by solving (36a) and (36b) under the same boundary conditions (equations (18)–(21b)) provided that the fictitious load distributions p
b
j
i
(
x
i
,
t
)
(
j
=
1
,
2,3
,
4
) are first established. These distributions can be determined following the procedure presented in [49] and employing the constant element assumption for the load distributions p
b
j
i along the L internal beam elements (as the numerical implementation becomes very simple and the obtained results are of high accuracy). Thus, the integral representations of the displacement components u
b
j
i (j
=
1
,
2,3
,
4) and their derivatives with respect to x
i when applied to the beam ends (0
,
l
i), together with the boundary conditions (18)–(21b), are employed to express the unknown boundary quantities u
b
j
i
(
ζ
i
), u
b
j
,
x
i
(
ζ
i
), u
b
j
,
x
x
i
(
ζ
i
), and u
b
j
,
x
x
x
i
(
ζ
i
) (ζ
i
=
0
,
l
i) in terms of p
b
j
i (j
=
1
,
2,3
,
4). Thus, the following set of 28 nonlinear algebraic equations for the ith beam is obtained as
(37)
[
E
b
11
i
0
0
0
0
E
b
22
i
0
0
0
0
E
b
33
i
0
0
0
0
E
b
44
i
]
{
d
b
1
i
d
b
2
i
d
b
3
i
d
b
4
i
}
+
{
0
D
b
1
i
n
l
0
0
D
b
2
i
n
l
0
0
D
b
3
i
n
l
0
0
0
}
=
{
0
α
b
3
i
0
0
β
b
3
i
0
0
γ
b
3
i
0
0
δ
b
3
i
}
,
where E
b
11
i and E
b
22
i
,
E
b
33
i
,
E
b
44
i are known matrices of dimension 4
×
(
N
+
4
) and 8
×
(
N
+
8
), respectively; D
b
1
i
n
l is a 2
×
1 and D
b
j
i
n
l (j
=
2,3) are 4
×
1 column matrices containing the nonlinear terms included in the expressions of the boundary conditions (equations (18)–(20b)); α
b
3
i and β
b
3
i, γ
b
3
i, δ
b
3
i are 2
×
1 and 4
×
1 known column matrices, respectively, including the boundary values of the functions a
b
3
i and β
b
3
i
,
β
-
b
3
i
,
γ
b
3
i
,
γ
-
b
3
i
,
δ
b
3
i
,
δ
-
b
3
i of (18)–(21b). Finally(38a)
d
b
1
i
=
{
p
b
1
i
u
^
b
1
i
u
^
b
1
,
x
i
}
T
,
(38b)
d
b
j
i
=
{
p
b
j
i
u
^
b
j
i
u
^
b
j
,
x
i
u
^
b
j
,
x
x
i
u
^
b
j
,
x
x
x
i
}
T
,
(
j
=
2,3
,
4
)are generalized unknown vectors, where(39a)
u
^
b
j
i
=
{
u
b
j
i
(
0
,
t
)
u
b
j
i
(
l
i
,
t
)
}
T
,
(
j
=
1
,
2,3
,
4
)
,
(39b)
u
^
b
j
,
x
i
=
{
∂
u
b
j
i
(
0
,
t
)
∂
x
i
∂
u
b
j
i
(
l
i
,
t
)
∂
x
i
}
T
,
(
j
=
1
,
2,3
,
4
)
,
(39c)
u
^
b
j
,
x
x
i
=
{
∂
2
u
b
j
i
(
0
,
t
)
∂
x
i
2
∂
2
u
b
j
i
(
l
i
,
t
)
∂
x
i
2
}
T
,
(
j
=
2,3
,
4
)
,
(39d)
u
^
b
j
,
x
x
x
i
=
{
∂
3
u
b
j
i
(
0
,
t
)
∂
x
i
3
∂
3
u
b
j
i
(
l
i
,
t
)
∂
x
i
3
}
T
,
(
j
=
2,3
,
4
)are vectors including the two unknown boundary values of the respective boundary quantities and p
b
j
i
=
{
(
p
b
j
i
)
1
(
p
b
j
i
)
2
⋯
(
p
b
j
i
)
L
}
T
(
j
=
1
,
2,3
,
4
) are vectors including the L unknown nodal values of the fictitious loads.Discretization of the integral representations of the unknown quantitiesu
b
j
i (j
=
1
,
2,3
,
4) and those of their derivatives with respect to x
i inside the beam x
i
∈
(
0
,
l
i
) and application to the L collocation nodal points yields(40a)
u
b
1
i
=
B
b
1
i
d
b
1
i
,
(40b)
u
b
1
,
x
i
=
B
b
1
,
x
i
d
b
1
i
,
(40c)
u
b
j
i
=
B
b
j
i
d
b
j
i
,
(
j
=
2,3
,
4
)
,
(40d)
u
b
j
,
x
i
=
B
b
j
,
x
i
d
b
j
i
,
(
j
=
2,3
,
4
)
,
(40e)
u
b
j
,
x
x
i
=
B
b
j
,
x
x
i
d
b
j
i
,
(
j
=
2,3
,
4
)
,
(40f)
u
b
j
,
x
x
x
i
=
B
b
j
,
x
x
x
i
d
b
j
i
,
(
j
=
2,3
,
4
)
,where B
b
1
i
,
B
b
1
,
x
i are L
×
(
L
+
4
) and B
b
j
i
,
B
b
j
,
x
i
,
B
b
j
,
x
x
i
,
B
b
j
,
x
x
x
i (j
=
2,3
,
4) are L
×
(
L
+
8
) known coefficient matrices, respectively.Applying (17a)–(17d) to the L collocation points and employing (40a)–(40f), 4
×
L nonlinear algebraic equations for each (ith) beam are formulated as(41a)
-
E
b
i
A
b
i
[
p
b
1
i
+
(
B
b
2
,
x
i
d
b
2
i
)
d
g
B
b
2
,
x
x
i
d
b
2
i
hhhhhhh
+
(
B
b
3
,
x
i
d
b
3
i
)
d
g
B
b
3
,
x
x
i
d
b
3
i
]
+
ρ
b
i
Α
b
i
T
1
q
¨
1
=
q
x
1
i
+
q
x
2
i
,
(41b)
E
b
i
I
z
i
p
b
2
i
-
(
N
b
i
)
d
g
B
b
2
,
x
x
i
d
b
2
i
-
ρ
b
i
I
z
i
B
b
2
,
x
x
i
d
¨
b
2
i
+
ρ
b
i
Α
b
i
B
b
2
i
d
¨
b
2
i
-
ρ
b
i
Α
b
i
(
B
b
1
i
d
¨
b
1
i
)
d
g
B
b
2
,
x
i
d
b
2
i
=
q
y
1
i
+
q
y
2
i
-
(
B
b
2
,
x
i
d
b
2
i
)
d
g
(
q
x
1
i
+
q
x
2
i
)
-
X
b
y
i
(
q
x
1
i
+
q
x
2
i
)
,
(41c)
E
b
i
I
y
i
p
b
3
i
-
(
N
b
i
)
d
g
B
b
3
,
x
x
i
d
b
3
i
-
ρ
b
i
I
y
i
B
b
3
,
x
x
i
d
¨
b
3
i
+
ρ
b
i
Α
b
i
B
b
3
i
d
¨
b
3
i
-
ρ
b
i
Α
b
i
(
B
b
1
i
d
¨
b
1
i
)
d
g
B
b
3
,
x
i
d
b
3
i
=
q
z
1
i
+
q
z
2
i
-
(
B
b
3
,
x
i
d
b
3
i
)
d
g
(
q
x
1
i
+
q
x
2
i
)
+
X
b
z
i
(
q
x
1
i
+
q
x
2
i
)
,
(41d)
E
b
i
C
S
i
p
b
4
i
-
G
b
i
I
t
i
B
b
4
,
x
x
i
d
b
4
i
+
ρ
b
i
I
p
i
B
b
4
i
d
¨
b
4
i
-
ρ
b
i
C
S
i
B
b
4
,
x
x
i
d
¨
b
4
i
=
e
y
1
i
q
z
1
i
+
e
y
2
i
q
z
2
i
-
e
z
1
i
q
y
1
i
-
e
z
2
i
q
y
2
i
+
Χ
b
w
i
(
q
x
1
i
+
q
x
2
i
)
,where (
N
b
i
)
d
g is a diagonal L
×
L matrix including the values of the axial forces of the ith beam; the symbol (
·
)
d
g indicates a diagonal L
×
L matrix with the elements of the included column matrix. The matrices X
b
z
i, X
b
y
i, X
b
w
i result after approximating the derivatives of m
b
y
j
i, m
b
z
j
i, m
b
w
j
i using appropriately central, backward, or forward differences. Their dimensions are also L
×
L. Moreover e
y
1
i, e
y
2
i, e
z
1
i, e
z
2
i are diagonal L
×
L matrices including the values of the eccentricities e
y
j
i, e
z
j
i of the components q
z
j
i, q
y
j
i with respect to the ith beam shear center axis (coinciding with its centroid), respectively.Employing (34a)–(34f) and (40a)–(40f), the discretized counterpart of (26a)-(26b), (27a)-(27b) and (28a)-(28b) at the L nodal points of each interface, is written as(42a)
Y
1
B
p
1
d
p
1
-
B
b
1
i
d
b
1
i
=
h
p
2
Y
1
B
p
3
,
x
d
p
3
+
h
b
i
2
B
b
3
,
x
i
d
b
3
+
b
f
i
4
B
b
2
,
x
i
d
b
2
+
(
φ
S
P
i
)
f
1
B
b
4
,
x
i
d
b
4
+
K
x
1
i
[
1
-
1
2
(
B
b
3
,
x
i
d
b
3
)
d
g
(
B
b
3
,
x
i
d
b
3
)
d
g
]
(
q
x
1
i
+
q
x
2
i
)
,
(42b)
Y
2
B
p
1
d
p
1
-
B
b
1
i
d
b
1
i
=
h
p
2
Y
2
B
p
3
,
x
d
p
3
+
h
b
i
2
B
b
3
,
x
i
d
b
3
-
b
f
i
4
B
b
2
,
x
i
d
b
2
+
(
φ
S
P
i
)
f
2
B
b
4
,
x
i
d
b
4
+
K
x
2
i
[
1
-
1
2
(
B
b
3
,
x
i
d
b
3
)
d
g
(
B
b
3
,
x
i
d
b
3
)
d
g
]
(
q
x
1
i
+
q
x
2
i
)
,
(43a)
Y
1
B
p
2
d
p
2
-
B
b
2
i
d
b
2
i
=
h
p
2
Y
1
B
p
3
,
y
d
p
3
+
h
b
i
2
B
b
4
i
d
b
4
i
+
K
y
1
i
[
1
-
1
2
(
B
b
4
,
x
i
d
b
4
)
d
g
(
B
b
4
,
x
i
d
b
4
)
d
g
]
(
q
x
1
i
+
q
x
2
i
)
,
(43b)
Y
2
B
p
2
d
p
2
-
B
b
2
i
d
b
2
i
=
h
p
2
Y
2
B
p
3
,
y
d
p
3
+
h
b
i
2
B
b
4
i
d
b
4
i
+
K
y
2
i
[
1
-
1
2
(
B
b
4
,
x
i
d
b
4
)
d
g
(
B
b
4
,
x
i
d
b
4
)
d
g
]
(
q
x
1
i
+
q
x
2
i
)
,
(44a)
Y
1
B
p
3
d
p
3
-
B
b
3
i
d
b
3
i
=
-
b
f
i
4
B
b
4
i
d
b
4
i
,
(44b)
Y
2
B
p
3
d
p
3
-
B
b
3
i
d
b
3
i
=
b
f
i
4
B
b
4
i
d
b
4
i
.Equations (30), (35a)–(35c), (37), and (41a)–(41d) together with continuity conditions (42a), (42b), (43a), (43b), (44a), and (44b) constitute a nonlinear system of algebraic equations with respect to q
x
1, q
x
2, q
y
1, q
y
2, q
z
1, and q
z
2 (interface forces) and d
p
i (i
=
1,2
,
3) and d
b
j
i (i
=
1
⋯
I), (j
=
1,2
,
3,4) (generalized unknown vectors of the plate and the beams). This system is solved using iterative numerical methods [52]. It is worth noting here that Y
1, Y
2 are position L
×
M matrices which convert the matrices B
p
i (i
=
1,2
,
3), B
p
3
,
x, and B
p
3
,
y into corresponding ones with dimensions L
×
M, appropriately referring to the nodal points of the two interface lines f
j
=
1
i, f
j
=
2
i, respectively. Moreover, K
x
1
i, K
x
2
i, K
y
1
i, and K
y
2
i are diagonal L
×
L matrices including the value of the flexibilities 1
/
k
x
j
i and 1
/
k
y
j
i
(
j
=
1,2
) of the arbitrary distributed shear connectors along the x
i and y
i directions, respectively.Finally, it is worth noting that beams placed along the boundary of the plate are treated as every other stiffening beam, since the lines of action of the integrated interface force componentsq
x
j
i, q
y
j
i, and q
z
j
i
(
j
=
1,2
) will also be internal ones, taking special care during the numerical evaluation of the line integrals in order to avoid their “near singular integral behaviour.” According to this, boundary elements that are very close to each other (distance smaller than their length) are divided in subelements, in each of which Gauss integration is applied [51].


## 4. Numerical Examples

On the basis of the analytical and numerical procedures presented in the previous sections, a FORTRAN program has been written and representative examples have been studied to demonstrate the effectiveness, wherever possible, the accuracy, and the range of applications of the proposed method. It is noted that the term “linear analysis” appearing in all of the following sections refers to the solution of the previously obtained system of equations neglecting all of the nonlinear terms.Example 1.
A rectangular plate (h
p
=
0.2
cm, E
p
=
E
b
=
3
×
1
0
7
kN
/
m
2, ρ
p
=
ρ
b
=
2.5
tn
/
m
3, ν
p
=
ν
b
=
0.2) with dimensions L
x
×
L
y
=
18
cm
×
9
cm subjected to a suddenly applied uniform load g
=
50
kN
/
m
2 and stiffened by a rectangular beam of 0.5
cm height and 1.0
cm width, eccentrically placed with respect to the plate center line (Figure 4), has been studied. The plate is clamped along its small edges, while the rest of the edges are free according to the transverse and inplane boundary conditions.
In Table1, the torsion I
t and warping C
S constants of the beam cross-section and the values of the primary warping function (
φ
S
P
)
j
(
j
=
1,2
) at the nodes of the two interface lines are presented.
The connection between the slab and the beam is accomplished using a linear distribution of shear connectors along each interface. The adopted relationship for the shear connectors’ stiffness is given as(45)
k
x
j
i
=
k
y
j
i
=
250
K
|
x
i
-
l
b
x
i
2
|
kN
/
m
2
(
j
=
1,2
)
,
where K is a dimensionless magnification factor. In Table 2 the obtained maximum deflections w
p of the first cycle of motion of the stiffened plate, at the middle of the free edge A, are shown for various values of the factor K, performing either a linear or a nonlinear analysis. In Figure 5 the time histories of the deflection w
p
(
t
) at the middle point A of the free edge of the examined stiffened plate, for the case of full connection between the plate and the beam, taking into account or ignoring geometrical nonlinearities, are presented as compared with those obtained from FEM solutions employing 8-noded hexahedral solid finite elements [25].
Moreover, in Figure6 the contour lines of the deflection w
p of the stiffened plate (full connection) at the time of maximum transverse displacement are presented as compared with those obtained from the aforementioned FEM solutions, ignoring (Figures 6(a) and 6(b)) or taking into account (Figures 6(c) and 6(d)) geometrical nonlinearities. In order to demonstrate the influence of the shear connectors in the dynamic behaviour of the stiffened plate, in Figure 7, the time histories of the deflection w
p
(
t
) at the same point A are presented, for various values of the factor K, performing either a linear or a nonlinear analysis. In Figure 8 the contour lines of the deflections w
p of the stiffened plate at the time of maximum displacement are presented, for various cases of connectors’ stiffness, ignoring (Figures 8(a), 8(c), and 8(e)) or taking into account (Figures 8(b), 8(d), and 8(f)) geometrical nonlinearities. Moreover, in Figure 9, the contour lines of the displacements u
p, v
p of the stiffened plate at the time of maximum deflection w
p are presented employing nonlinear dynamic analysis.Table 1
Torsion, warping constants, and primary warping function at the interface lines of the stiffening beam of Example1.
h
b (cm)
I
t (cm4)
C
S (cm6)
(
ϕ
S
P
)
f
1 (cm2)
(
ϕ
S
P
)
f
2 (cm2)
0.5
2.859
E
-
02
3.175
E
-
04
-
4.979
E
-
02
4.979
E
-
02Table 2
Maximum deflectionw
p (cm) at point A of the first cycle of motion of the stiffened plate for various values of K of Example 1.
K
Linear analysis
Nonlinear analysis
Full connection
0.484
0.391
100.0
0.488
0.395
10.0
0.495
0.401
1.0
0.521
0.419
0.1
0.575
0.449
0.01
0.696
0.494
0.005
0.725
0.501Plan view (a) and section a-a (b) of the stiffened plate of Example1.
(a)
(b)Figure 5
Time history of deflectionw (cm) at the middle point A of the free edge of the plate of Example 1.Contour lines ofw
p of the stiffened plate of Example 1, employing the present study (a, c) and a FEM [25] solution using solid elements (b, d) at the time of maximum transverse displacement.
(a)
w
p
max
⁡
=
0.484
cm at t
=
1.19 msec of linear AEM analysis
(b)
w
p
max
⁡
=
0.491
cm at t
=
1.20 msec of linear FEM [25] analysis
(c)
w
p
max
⁡
=
0.394
cm at t
=
1.08 msec of nonlinear AEM analysis
(d)
w
p
max
⁡
=
0.387
cm at t
=
1.06 msec of nonlinear FEM [25] analysisFigure 7
Time history of the deflectionw (cm) at the middle point A of the free edge of the stiffened plate of Example 1 for various values of connectors’ stiffness.Contour lines ofw
p of the stiffened plate of Example 1 at the time of maximum transverse deflection ignoring (a, c, e) or taking into account geometrical nonlinearities (b, d, f).
(a)
w
p
max
⁡
=
0.725 cm at t
=
1.47 msec of linear AEM analysis for K
=
0.005
(b)
w
p
max
⁡
=
0.501 cm at t
=
1.23 msec of nonlinear AEM analysis for K
=
0.005
(c)
w
p
max
⁡
=
0.521 cm at t
=
1.22 msec of linear AEM analysis for K
=
1.0
(d)
w
p
max
⁡
=
0.419 cm at t
=
1.11 msec of nonlinear AEM analysis for K
=
1.0
(e)
w
p
max
⁡
=
0.495 cm at t
=
1.20 msec of linear AEM analysis for K
=
10.0
(f)
w
p
max
⁡
=
0.401 cm at t
=
1.08 msec of nonlinear AEM analysis for K
=
10.0Contour lines ofu
p, v
p of the stiffened plate of Example 1 at the time of maximum transverse deflection w
p taking into account geometrical nonlinearities.
(a)
u
p
max
⁡
=
4.5
×
10
-
3 cm at t
=
1.23 msec for K
=
0.005
(b)
v
p
max
⁡
=
9.3
×
10
-
3 cm at t
=
1.23 msec for K
=
0.005
(c)
u
p
max
⁡
=
4.3
×
10
-
3 cm at t
=
1.11 msec for K
=
1.0
(d)
v
p
max
⁡
=
6.4
×
10
-
3 cm at t
=
1.11 msec for K
=
1.0
(e)
u
p
max
⁡
=
4.0
×
10
-
3 cm at t
=
1.08 msec for K
=
10.0
(f)
v
p
max
⁡
=
6.0
×
10
-
3 cm at t
=
1.08 msec for K
=
10.0Example 2.
The linear and nonlinear response of a rectangular plate having a central stiffener, as shown in Figure10, has been studied (E
=
68.9
GPa, v
=
0.30, ρ
=
2.67
tn
/
m
3). The plate is subjected to a suddenly applied uniform load of g
=
300
kPa, while its boundaries are simply supported with restraint against inplane motion. In Figure 11 the linear and nonlinear time history response of the deflection of the half panel center A is depicted. In this figure the obtained results are compared with those presented by Jiang and Olson [53] employing conventional finite strip method, by Koko [54] employing super finite element method, and by Sheikh and Mukhopadhyay [22] employing the spline finite strip method.
As it can easily be observed, there is significant agreement between the results concerning the linear dynamic analysis while deviation between the different types of analysis is observed in nonlinear dynamic analysis, which has been already mentioned by the previous investigators [22, 53, 54].Plan view (a) and section a-a (b) of the stiffened plate of Example2.
(a)
(b)Time history of linear (a) and nonlinear (b) response of the deflectionw (mm) of the half panel center A of the stiffened plate of Example 2.
(a)
(b)Example 3.
A rectangular plate with dimensionsL
x
×
L
y
=
1.8
m
×
0.9
m (Figure 12) stiffened by an I cross-section beam (Figure 13) eccentrically placed with respect to the plate centre line, clamped along all its edges, and subjected to a suddenly applied uniform load g
=
750
kN
/
m
2 has been studied (h
p
=
0.02
m, E
p
=
E
b
=
68.9
×
1
0
6
kN
/
m
2, ρ
p
=
ρ
b
=
2.67
tn
/
m
3, ν
p
=
ν
b
=
0.3). In Table 3, the torsion I
t and warping C
S constants of the beam cross-section and the values of the primary warping function (
φ
S
P
)
j
(
j
=
1,2
) at the nodes of the two interface lines are presented. In Figure 14 the time histories of the maximum deflection w
p
(
t
) of the plate for the cases of linear and nonlinear dynamic analysis are presented as compared with those obtained from FEM solutions using 8-noded hexahedral solid finite elements [25]. Moreover, in Figure 15, the contour lines of the displacement w
p of the stiffened plate at the time of maximum transverse displacement are presented as compared with those obtained from the aforementioned FEM solutions, ignoring (Figures 15(a) and 15(b)) or taking into account (Figures 15(c) and 15(d)) geometrical nonlinearities. The convergence of the results between the two methods is noteworthy. Finally, in Figure 16, the contour lines of the corresponded moments M
x, M
y at the time of maximum transverse displacement are presented.Table 3
Torsion, warping constants, and primary warping function at the interface lines of the stiffening beam.
I
t (m4)
C
S (m6)
(
ϕ
S
P
)
f
1 (m2)
(
ϕ
S
P
)
f
2 (m2)
6
.
984
E
-
08
3.374
E
-
09
-
9.94
E
-
04
9.94
E
-
04Plan view (a) and section a-a (b) of the stiffened plate of Example3.
(a)
(b)Figure 13
Cross-section of the stiffener of Example3.Figure 14
Time history of the maximum displacementw (mm) of the stiffened plate of Example 3.Contour lines of deflectionw
p of the stiffened plate of Example 3, employing the present study (a, c) and a FEM [25] solution using solid elements (b, d) at the time of maximum transverse deflection.
(a)
w
p
max
⁡
=
24.5 mm at t
=
2.59 msec of AEM linear analysis
(b)
w
p
max
⁡
=
25.8 mm at t
=
2.65 msec of FEM linear analysis
(c)
w
p
max
⁡
=
23.0 mm at t
=
2.54 msec of AEM nonlinear analysis
(d)
w
p
max
⁡
=
22.6 mm at t
=
2.40 msec of FEM nonlinear analysisContour lines of momentsM
x and M
y (kNm/m) at the time of maximum transverse displacement.
(a)
M
x at t
=
2.59 msec of AEM linear analysis
(b)
M
x at t
=
2.54 msec of AEM linear analysis
(c)
M
y at t
=
2.59 msec of AEM nonlinear analysis
(d)
M
y at t
=
2.54 msec of AEM nonlinear analysisTaking into account the aforementioned convergence of the results between the two methods (AEM, FEM), the importance of the reduction of calculation time when employing the proposed method should be highlighted. More specifically, for the determination of the beam’s response to the dynamic loading, a personal computer of 8 Gb RAM and Intel I7 processor of 4 cores with maximum clock speed equal to 3.8 GHz was used. The time step used for the nonlinear analysis is 1 microsecond (μsec) and for a full circle response the calculation time employing the proposed method was 20 minutes as compared with FEM analysis which lasted 50 minutes.


## 5. Concluding Remarks

A general solution for the geometrically nonlinear dynamic analysis of plates stiffened by arbitrarily placed parallel beams of arbitrary doubly symmetric cross-section subjected to arbitrary loading is presented. The proposed model takes into account the nonuniform distribution of the interface shear forces and the nonuniform torsional response of the beams. The main conclusions that can be drawn from this investigation are as follows.(a)
The proposed model permits the dynamic response of stiffened plates subjected to arbitrary loading, while both the number and the placement of the parallel beams are also arbitrary (eccentric beams are included). The plate and the beams are supported by the most general boundary conditions including elastic support or restraint.
(b)
The adopted model permits the evaluation of the longitudinal and transverse inplane shear forces at the interfaces between the plate and the beams in the geometrically nonlinear analysis of the stiffened plate, the knowledge of which is very important in the design of shear connectors in stiffened structures.
(c)
The nonuniform torsion in which the stiffening beams are subjected is taken into account by solving the corresponding problem and by comprehending the arising twisting and warping in the corresponding displacement continuity conditions. The distributed warping moment arising from the nonuniform distribution of longitudinal inplane forces is also taken into account.
(d)
The accuracy of the results and the validity of the proposed model are noteworthy.
(e)
The influence of geometrical nonlinearity on the deformation of the examined stiffened plates is remarkable. In the case of immovable inplane boundary conditions, the displacements decrement can be verified.
(f)
The increment of the deflection with the decrement of the connectors’ stiffness is easily verified, while this decrement results in a more pronounced influence of the geometrical nonlinearity on the response of the stiffened plate.
(g)
The developed procedure retains most of the advantages of a BEM solution over a FEM approach, although it requires domain discretization.



---

*Source: 942763-2014-09-14.xml*