#### 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.

#### 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 elasticity , volume mass density , shear modulus , and Poisson’s ratio having constant thickness and occupying the two-dimensional multiple connected region of the plane bounded by the piecewise boundary curves. The plate is stiffened by a set of arbitrarily placed parallel beams of arbitrary doubly symmetric cross-section of area and length . The material of the beams is considered to be homogeneous, isotropic, and linearly elastic with modulus of elasticity , mass density , shear modulus , and Poisson’s ratio . Therefore, material nonlinearities are not considered in the following analysis.

For the sake of convenience the axis is taken parallel to the beams of length which may have either internal or boundary point supports. The stiffened plate is subjected to the arbitrary lateral dynamic load , , . 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 for the analysis of the plate and local coordinate ones corresponding to the centroid axes of each beam are employed (Figure 1).

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 th beam results in line forces per unit length in all directions in two interface lines, which are denoted by , , and (), encountering in this way the nonuniform distribution of the interface transverse shear forces , which in previous models [43] was ignored. The aforementioned integrated tractions result in the loading of the th 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.

**(a)**

**(b)**

*(a) In the Plate (at the Traces of the Two Interface Lines ** of the **th Plate-Beam Interface)*(i)An inplane line body force at the middle surface of the plate.(ii)An inplane line body force at the middle surface of the plate.(iii)A lateral line load .(iv)A lateral line load due to the eccentricity of the component from the middle surface of the plate. is the bending moment.(v)A lateral line load due to the eccentricity of the component from the middle surface of the plate. is the bending moment.

*(b) In Each (**th) Beam (** System of Axes)*(i)An axially distributed line load along the beam centroid axis .(ii)A transversely distributed line load along the beam centroid axis .(iii)A perpendicularly distributed line load along the beam centroid axis .(iv)A distributed bending moment along local beam centroid axis due to the eccentricities of the components from the beam centroid axis. are the eccentricities.(v)A distributed bending moment along local beam centroid axis due to the eccentricities of the components from the beam centroid axis. and are the eccentricities.(vi)A distributed twisting moment along local beam shear center axis due to the eccentricities and of the components and from the beam shear center axis, respectively. and and are the eccentricities.(vii)A distributed warping moment ,along local beam shear center axis, was ignored in previous models [39, 43]. 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 th interface line of the th 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 aswhere , , , , are the time dependent inplane and transverse displacement components of an arbitrary point of the plate and , , and ,, 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 asSubstituting (1a)–(1c) and (2a)–(2d) to the stress-strain relations defined by the Hooke’s law
the nonvanishing components of the second Piola-Kirchhoff stress tensor are obtained asSubsequently, integrating the stress components over the plate thickness, the stress resultants acting on the plate are written aswhere and 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 aswhere is the Dirac’s delta function in the direction. Employing relations (5a)–(5f), the governing differential equations (6a)–(6c) in the domain can be expressed in terms of the displacement components asThe governing differential equations (7a)–(7c) are also subjected to the pertinent boundary conditions of the problem at handand to the initial conditionswhere , , , and are functions specified at the boundary ; are functions specified at the corners of the plate; , , and are the initial deflection and velocity of the points of the middle surface of the plate; , and , are the boundary membrane displacements and forces in the normal and tangential directions to the boundary, respectively; and 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 to the boundary and the arc length ) are written asin which is the curvature of the boundary. Finally, is the discontinuity jump of the twisting moment at the corner of the plate, while along the boundary is given by the following relation:
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 corners* k*. It is apparent that all types of the conventional boundary conditions can be derived from these equations by specifying appropriately the functions , , , and () (e.g., for a clamped edge it is , ).

*(b) For Each (**th) Beam.* Each beam undergoes transverse deflection with respect to and axes and axial deformation and nonuniform angle of twist along 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 aswhere , , and are the axial and transverse displacement components with respect to the system of axes; , , and are the corresponding components of the centroid ; and are the angles of rotation of the cross-section due to bending, with respect to its centroid; denotes the rate of change of the angle of twist regarded as the torsional curvature and 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 asEmploying 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 aswhere is the primary twisting moment [15] resulting from the primary shear stress distribution; is the warping moment due to torsional curvature. Furthermore and are the principal moments of inertia; is the polar moment of inertia, while and are the torsion and warping constants of the th beam with respect to the cross-section’s shear center (coinciding with its centroid), respectively, given as [46] On the basis of Hamilton’s principle, the differential equations of motion in terms of displacements are obtained asSubstituting the expressions of the stress resultants of (14a)–(14e) in (16a)–(16d), the differential equations of motion are obtained asMoreover, the corresponding boundary conditions of the th beam at its ends are given as
and the initial conditions aswhere the angles of rotation of the cross-section due to bending , are given from (12d) and (12e), , and , are the reactions and bending moments with respect to , axes, respectively, which after applying the aforementioned simplifications are given asand , are the torsional and warping moments at the boundaries of the beam, respectively, given asFinally, () are functions specified at the th beam ends (). 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 , , , , , , , , , , , and . Six additional equations are required, which result from the displacement continuity conditions in the directions of , , and local axes along the two interface lines of each (th) 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 of local axis: In the direction of local axis: In the direction of local axis: where 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 th interface line of the th plate-beam interface and , are the stiffness of the arbitrarily distributed shear connectors along the and directions, respectively. It is noted that and can represent any linear or nonlinear relationship between the inplane interface forces and the interface slip in the corresponding direction. In all of the aforementioned equations the values of the primary warping function 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 Components , ,

According to the precedent analysis, the large deflection analysis of the plate becomes equivalent to establishing the inplane displacement components and having continuous partial derivatives up to the second order with respect to , , and the deflection having continuous partial derivatives up to the fourth order with respect to , , and all displacement components having derivatives up to the second order with respect to . 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 , , and applying the Laplacian operator to , and the biharmonic operator to yieldsEquations (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 are first established. These distributions can be determined using BEM. Following the procedure presented in [47], the unknown boundary quantities , and (, ) can be expressed in terms of () after applying the integral representations of the displacement components and their derivatives with respect to to the boundary of the plate.

By discretizing the boundary of the plate into boundary elements and the domain into 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 and of the boundary conditions (8a)–(8e) to the boundary nodal points results in a set of nonlinear algebraic equations relating the unknown boundary quantities with the fictitious load distributions () that can be written as where to are rectangular known matrices including the values of the functions () of (8a)–(8d); , , , and are known column matrices including the boundary values of the functions , of (8a)–(8d); (), , and () are rectangular known coefficient matrices resulting from the values of kernels at the boundary elements of the plate; , , and are rectangular known matrices originating from the integration of kernels on the domain cells of the plate; () are and is 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 , appearing in (8a)–(8d), are approximated employing appropriate central, backward, or forward finite difference schemes. Finallyare generalized unknown vectors, where are vectors including the unknown boundary values of the respective boundary quantities and , (), are vectors containing the unknown nodal values of the fictitious loads at the domain cells of the plate. In the case that the boundary has free or transversely elastically restrained corners, additional equations must be satisfied together with (30). These additional equations result from the application of the corner condition (8e) on the corners following the procedure presented in [48].

Discretization of the integral representations of the displacement components and their derivatives with respect to , [49] and application to the domain nodal points yieldswhere , , , , , () are and , , , , , are 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 nodal points inside , yielding the following equations:where , , and are vectors of dimension 2 including the unknown , , () interface forces; is the total number of the nodal points at each plate-beam interface; is a position matrix which converts the vectors , , into corresponding ones with length ; the symbol indicates a diagonal matrix with the elements of the included column matrix. Matrices and of dimension result after approximating the bending moment derivatives of , , respectively, using appropriately central, backward, or forward differences.

##### 3.2. For the Beam Displacement Components , , and and for the Angle of Twist

According to the precedent analysis, the nonlinear dynamic analysis of the th beam reduces in establishing the displacement components and , , having continuous derivatives up to the second order and up to the fourth order with respect to , respectively, and also having derivatives up to the second order with respect to . 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).

Let , , , and be the sought solutions of the problem represented by (17a)–(17d). Differentiating with respect to these functions, two and four times, respectively, yieldsEquations (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 are first established. These distributions can be determined following the procedure presented in [49] and employing the constant element assumption for the load distributions along the 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 () and their derivatives with respect to when applied to the beam ends (), together with the boundary conditions (18)–(21b), are employed to express the unknown boundary quantities , , , and () in terms of (). Thus, the following set of 28 nonlinear algebraic equations for the th beam is obtained as where and are known matrices of dimension and , respectively; is a and () are column matrices containing the nonlinear terms included in the expressions of the boundary conditions (equations (18)–(20b)); and , , are and known column matrices, respectively, including the boundary values of the functions and of (18)–(21b). Finallyare generalized unknown vectors, whereare vectors including the two unknown boundary values of the respective boundary quantities and are vectors including the unknown nodal values of the fictitious loads.

Discretization of the integral representations of the unknown quantities () and those of their derivatives with respect to inside the beam and application to the collocation nodal points yieldswhere are and