Extension of the Wittrick-Williams Algorithm for Free Vibration Analysis of Hybrid Dynamic Stiffness Models Connecting Line and Point Nodes

: This paper extends the Wittrick-Williams (W-W) algorithm for hybrid dynamic stiffness (DS) models connecting any combinations of line and point nodes. The principal novelties lie in the development of both the DS formulation and the solution technique in a sufﬁciently systematic and general manner. The parent structure is considered to be in the form of two dimensional DS elements with line nodes, which can be connected to rigid/spring point supports/connections, rod/beam point supports/connections, and point connections to substructures. This is achieved by proposing a direct constrain method in a strong form which makes the modeling process straightforward. For the solution technique, the W-W algorithm is extended for all of the above hybrid DS models. No matrix inversion is needed in the proposed extension, making the algorithm numerically stable, especially for complex built-up structures. A mathematical proof is provided for the extended W-W algorithm. The proposed DS formulation and the extended W-W algorithm are validated by the FE results computed by ANSYS. This work signiﬁcantly extends the application scope of the DS formulation and the W-W algorithm in a methodical and reliable manner, providing a powerful eigenvalue analysis tool for beam-plate built-up structures.


Introduction
Built-up structures in engineering are often composed of plates, beams and other elements through a variety of imposed constraints or connections [1][2][3]. The applicability of the built-up structures includes, but is not limited to, solar panels, building floors, ships, aircraft, armored vehicles, cars, robots, and circuit boards, among others. Vibration analysis is an essential requirement for the design of these built-up structures. Moreover, supports/connections can significantly affect the dynamic behavior of built-up structures, such as natural frequency, transmission path, and vibro-acoustic characteristics. If these constraints or connections are not properly considered in the design, inaccurate dynamic analysis may result in undesirable phenomena such as resonance or structural fatigue. The supports/connections of beam-plate combination of built-up structures can be divided into continuous (linear) supports/connections and discrete (point) supports/connections. Point supports/connections are very common in built-up structures in engineering. Figure 1 describes a range of common point supports/connections, which can be divided into two main categories: (a) Point supports (PS), including rigid support, spring support, rod support and beam support, see Figure 1a, and (b) point connections (PC), including rigid connection, spring connection, rod connection and beam connection, see Figure 1b. In the modeling of built-up structures shown in Figure 1, various supports and connections for plate structures can be generally idealized as points, springs, beams and rods connected to the plate. That is, the line nodes are connected with various forms of point nodes (PS/PC) in the built-up structure.
There are two main types of modeling methods for the built-up structures under consideration, namely, numerical methods and analytical/semi-analytical methods. Numerical methods include finite element method (FEM) [4][5][6][7][8], finite difference method [9][10][11], modal constraint method [12][13][14], Galerkin method [15,16], etc. However, numerical methods often fail to achieve high accuracy and computational efficiency at the same time. For example, Mello et al. [4] developed a finite element model for built-up systems made of concrete slabs and steel beams and carried out a dynamic analysis and showed that accuracy of result and computational efficiency run encounter. Machado et al. [5] developed a dynamic finite element formulation (plate, bar and interface) to model building floors, and evaluated the influence of connection stiffness on the dynamic properties. Altabey [6] used a finite strip transition matrix for the response analysis of laminated structures with four classical boundary conditions under different elastic constraint coefficients. Recently, Friswella and Wang [7,8] used FEM to model plates with additional supports/mass, and investigated the optimal position and minimum constraint stiffness of flexible PS to attenuate the natural frequencies of the plate structures, but their modeling was quite complicated. By contrast, analytical or semi-analytical modeling methods can be both computationally accurate and efficient at the same time, for example, the Ritz method [17][18][19][20][21][22][23][24], finite strip method [25][26][27][28][29], superposition method [16,[30][31][32][33][34][35], and Lagrange multiplier method [36][37][38]. In the analytical or semi-analytical method, modeling supports/connections of beamplate built-up structures can also be divided into continuous (line supports/connections) and point supports/connections like the numerical methods. For instance, continuous (line) supports/connections (including beam stiffeners) are essentially modeled as direct supports or connections to the line nodes. It should be noted that for uniformly distributed elastic supports/rigid supports, the Ritz method [17][18][19][20][21], Fourier series based analytical method [39,40], and finite strip method [28,29] have routinely been used. However, the modeling of arbitrarily spaced non-uniform continuous support are no-doubt difficult [41][42][43][44][45][46][47]. Similarly, the beam-stiffened plate that is generally modeled as line nodes can be connected to the line nodes of plate boundaries based on different methods, such as generalized differential quadrature method [48], superposition method [31], mode-based approach [49], Ritz method based on modified Fourier series [50], and double finite cosine integral transformation method [51].
However, the analytical modeling for a built-up structure with a line node connected to point nodes, shown in Figure 1, is significantly more challenging than those of continuous support/connection modeling. There are few existing works in this area, which can be classified into two groups, which are explained in detail as follows: Group (i) analytical modeling of line nodes subjected to point supports. For instance, Gorman and Singal [32] proposed the analytical free vibration solutions of rectangular plates constrained by bolts or spot welding based on superposition method. Huang and Thambiratnam [27] proposed a finite strip element method combined with spring system for the free vibration analysis of elastic intermediate support plate. Dozio [20] extended the trigonometric Ritz method for the vibration characteristics of plate structures with midline and point supports. Tripathy and Suryanaraya [52,53] studied the vibration and buckling characteristics of welded rectangular plate by the flexible function method. Recently, Li et al. [33,34] proposed a symplectic superposition method for the free vibration of rectangular plate with multi-point supports. Group (ii) Analytical modeling of line nodes with point connected nodes to flexible beams. For instance, Chiba and Yoshida [22] used the Ritz method for free vibration analysis of a cantilever rectangular plate coupled to beam structures. Garusi and Tralli [54] proposed a new transition element based on the FEM and mixed stress method to model solid-to-beam and plate-to-beam connections. Cao et al. [2,3] used the global modal method combined with the Ritz method to analyze the vibration characteristics of spacecraft flexible joints, where the high-aspect-ratio plate was modeled as a beam, but the coupling of bending and torsional vibration was ignored. Furthermore, some electronic microstructure components, sometimes called the island-bridge model, are essentially plated like island connected to beam-like bridges, have received attention of many investigators [1,[55][56][57]. Nevertheless, majority of these investigations are focused on static and buckling analysis only. It can be seen from the above literature review that most analytical solutions are for continuous support (line nodes connected to line nodes), while analytical solutions for built-up structures with discrete connections (line nodes connected to point nodes) are relatively rare, and the modeling process is somehow cumbersome. Moreover, modeling for different support or connection often requires separate and different formulations, which are not really convenient for straightforward application in engineering structures.
Against the above background, the DSM [58] is a powerful alternative modeling method for built-up structures. This is because first, the DSM is an exact analytical method, in which exact individual elements are based on exact shape functions, and thus is capable of computing exact results by using very few degrees of freedom (DOF) within the whole frequency ranges. Secondly, there is a comprehensive group of different elements in DSM, and importantly, the elements can be assembled directly to model complex builtup structures. When the analytical DSM is applied to modal analysis, there are many eigenvalues solution techniques, such as the determinant method [59][60][61][62][63]. However, the determinant method is inefficient, and some natural frequencies may be missed. On the contrary, the Wittrick-Williams (W-W) algorithm [64] is an accurate, efficient and robust method to compute eigenvalues from analytical DS models, which ensures that no natural frequencies will be missed. The powerful combination of DS formulations and the W-W algorithm has been applied to a wide range of problems in the literature, such as the free vibration analyses of advanced rods [65], tapered beams [66], lattice microstructures [67][68][69], cracked frames [70][71][72], multi-body structures [73][74][75], plates [76][77][78][79][80], membranes [81] and recently for uncertain structures [82,83]. Once the DS elements have been developed, the solution of built-up structures containing 1D and 2D elements with the application of W-W algorithm is quite straightforward. This is because the 1D elements are assembled directly by point nodes, and the 2D elements are assembled directly via line nodes.
There are few reported research on DS formulation and the associated W-W algorithm devoted to built-up structures with line nodes connected to point nodes. In this respect, Williams and Anderson [84] proposed a DS formulation method and applied the Wittrick-Williams algorithm based on the Lagrange Multipliers. This approach was subsequently used for free vibration [38] and wave propagation [85] of plate elements with line nodes and with rigidly/elastically point supports. On the other hand, Lam et al. [37] investigated the buckling of two overlap-jointed plate elements connected by bolts and spot welding. Interestingly, Powell et al. [86] gave a different perspective by applying the method for efficient multi-level substructuring with constraint for buckling and vibration of prismatic plate assemblies. The related formulation and the W-W algorithm have also been extended to model hybrid models by connecting line nodes of DS model with the point nodes of FE model, see [87,88]. However, all of the above works are based on dynamic flexibility formulation rather than DS formulations. The calculation process involved cumbersome matrix inversion, which may become numerically instable for larger systems with more degrees of freedom. Moreover, the above research was by and large applied to plates with point supports, and the formulation and the associated W-W algorithm application were not extended to plate elements with line nodes connected to point nodes of DS elements of bars and beams. This is a significant omission in the literature.
Motivated by the work of Williams and Anderson [84], the main purpose of this paper is to propose a unified and universal DS modeling method that extends the W-W algorithm for built-up structures with line nodes connected to point nodes, including eight types of different supports/connections, as illustrated in Figure 1.
DS formulation for 1D elements, 2D elements and the classical W-W algorithm are not given detailed coverage here because they are available in the literature. The paper proceeds as follows. First, the DS formulation for 1D elements, 2D elements and the classical Wittrick-Williams algorithm are briefly reviewed in Section 2. Then, the new DS formulation with extended W-W algorithm for built-up structures with line nodes connected to point nodes are presented in Section 3. Next, the proposed formulation and the extended W-W algorithm are validated against the FEM results by modal analysis of several representative structures in Section 4. Finally, conclusions of this work are drawn in Section 5. This paper is basically a significant extension of the work of Williams and Anderson [84]. Additionally, the work in this paper can be applied in a much wider context of built-up structures, including the coupling effects of 1D DS elements and 2D elements in the general case. An important attribute of this paper is that the modeling procedure is more direct and convenient with much better numerical stability properties.

Dynamic Stiffness Formulation for 1D and 2D Elements and the Classical Wittrick-Williams Algorithm
In Section 2.1 below, the general procedure for the DS formulation for 1D and 2D elements is presented. Then, the DS matrix for 2D plate element, 1D bar and beam elements, and the massless spring element are presented, respectively, in Sections 2.2-2.4. Finally, the classical Wittrick-Williams algorithm is briefly discussed in Section 2.5.

General Procedure for the DS Formulation for 1D Element, 2D Element
The general procedure for the DS formulation for a structural element (either 1D or 2D) is provided as follows. Based on Hamilton's principle or other suitable methods, such as Newton's Laws and Lagrange's equation, the governing differential equations (GDEs) describing the structural motion in the time domain are obtained in the general form as follows L(u, t) = 0 where L is a linear differential operator, u is displacement vector , t is time. Then, based on harmonic oscillation assumption by letting u(x, t) = U(x)e iωt , where U(x) is displacement amplitude vector, ω stands for circular frequency, the GDEs can be written as where L 1 is the linear differential operator in the frequency domain. Subsequently, the exact general solution of Equation (2) can be derived U(x) = N(x)C, where N(x) is the exact shape functions and C is the unknown vector of constants. By eliminating constant vector C using boundary conditions of force and displacement, the elemental DS matrix can be obtained as where f e and d e are, respectively, the generalized amplitudes of boundary force and displacement of the structure element in the frequency domain, and K e (ω) is the DS matrix. Next, the DS formulation of 2D plate elements and 1D beam/spring elements are briefly reviewed as follows.

The DS Formulation for 2D Plate Element
The GDE of classical plate theory of a thin Levy-type plate element ( Figure 2) in the frequency domain is [76]  Considering the Levy-type boundary conditions, the general solution of Equation (4) can be written as The force and displacement boundary conditions along the nodal edges of the plate element of half-wavenumber m are, respectively: Finally, the DS matrix of plate can be written in the following general form [76,77]: i.e., where Q i and M i (i = 1, 2) represent generalized amplitudes of boundary forces, w and φ represent generalized amplitudes of boundary displacements, and the analytical expressions of k 11 ∼k 24 in Equation (7) are given in Appendix A, also see Equations (A1) and (A2).

The DS Formulation of 1D Rod and Beam Element
This section gives the DS formulation for a beam element based on the classical rod theory for axial vibration and the Euler-Bernoulli theory for transverse vibration, see Figures 3 and 4.
The GDE of a rod in the frequency domain can be written as follows [73]: is the amplitude of axial vibration and ξ = x/l. Now, referring to Figure 3, the boundary conditions for displacements and forces at both ends can be applied as follows: The DS matrix of the rod can be written in the following form: i.e., where the analytical expressions of the elements e 1 and e 2 in Equation (12) are given in Appendix B. The frequency domain GDE of a Euler-Bernoulli beam for bending vibration is given as follows [73]: where λ = mω 2 l 4 EI , W(ξ) is the amplitude of transverse deflection, ξ = x/l. Now, referring to Figure 4, the boundary conditions for displacements and forces at both ends of the beam can be applied as follows: The DS matrix of a beam under bending vibration can now be written in the following general form: i.e., T denotes generalized amplitudes of boundary forces and T denotes generalized amplitude of boundary displacements. F y and M represent the shear force and bending moment, respectively, at the two ends denoting point nodes; w represents vertical or bending displacement, and the θ represents the bending rotation. Note that the 4 × 4 stiffness matrix shown in Equation (17) must be symmetric, d 1 ∼d 6 are the DS coefficients for beam vibration, which are given in Appendix C.

Massless Spring
Let us assume that the ith and the jth point nodes of a structure are connected by a massless elastic spring with stiffness k p , as shown in Figure 5. Here, k p stands for either translational or rotational stiffness. The displacement equation for the ith and jth point nodes will then be where d i and d j are the generalized amplitudes of displacements for the ith and the jth point nodes, respectively; f a i and f a j are the forces acting onto the corresponding point nodes due to the point elastic constraint with stiffness k p . Therefore, Equation (19) can be rewritten in the following matrix form as

Classical Wittrick-Williams Algorithm in Original Form
Once the global dynamic stiffness matrix of the whole structure is formulated, the W-W algorithm can be applied to compute the natural frequencies ω of the structure. The algorithm is based on the computation of the mode count J when ω is lower than the trial where [K(ω * )] is the DS matrix of overall structure, and s{[K(ω * )]} is the number of negative diagonal elements after upper triangular transformation by using Gauss elimination of [K(ω * )]. J 0 is given by where J 0 (ω * ) is the number of natural frequencies between ω = 0 and ω = ω * when the nodes of all the elements are fully clamped. The expressions for J 0m of plate, rod and beam are, respectively, given by Equations (A3), (A5) and (A8) (see Appendices A-C).

The DS Formulation and Extended Wittrick-Williams Algorithm for Built-Up Structures with Line Nodes Connected to Point Nodes
This section presents a direct constraint method for the connection of the DS matrices of 1D and 2D elements (see Section 3.1), and provides the extended W-W algorithm with the corresponding formulation procedure of hybrid DS modeling (see Section 3.2). As it will be seen that this useful extension is far from trivial, and efforts expended are very considerable.

The Direct Constraint Method
The proposed DS formulations and the associated Wittrick-Williams algorithm are intended for built-up structures composed of either 1D elements connected via point nodes or 2D elements connected via line nodes. However, for built-up structures where 2D elements with line nodes connected to point nodes require new DS formulation and the associated Wittrick-Williams algorithm application. In this section, a direct constraint method is proposed to connect a line node of a plate structure to a point node in a sufficiently general manner covering all possible cases, see Figure 1. The procedure is essentially based on a strong formulation using the equilibrium conditions at the point connections.
Based on the premises of the DS method, the generalized displacement d(ξ) along the ith line DOF (CPT, ξ ∈ [0, L]) can be written as where Therefore, the generalized displacement d ip at point ξ = ξ p on the ith line DOF can be written as where e(ξ p ) is a vector obtained by evaluating e(ξ) at point ξ = ξ p . If there are P point connections on the ith line node to be considered, then where For all of the point connections on all of the line DOF of the plate assembly, the corresponding equation can be written in the form where When the plate assembly is under free vibration condition with the specified boundary conditions and point connections, a reciprocal relation can be obtained between the point forces and the boundary forces. To be more specific, if concentrated point generalized forces f ip (p = 1∼P) are acting on P points of the ith line DOF, then these point forces will cause additional generalized force vector f a i to the ith line DOF becomes The above equation can be obtained by integrating the Dirac Delta functions which represent the concentrated generalized forces f ip on the ith line DOF. Similarly, the corresponding equation for the whole plate assembly becomes The boundary force vector f becomes the summation of the plate response f Pl = K Pl d Pl and the additional boundary response resulting from these point connections f a .
Therefore, the equilibrium condition can be written as namely, where K plm represents the DS matrix of the plate [38,84], and the assembly of the plate element refers to [77], m stands for the half-wavenumber of a plate and I is the unit matrix.
It should be pointed out that the method proposed in this paper can be directly applied for composite plates [89] and plate assemblies with arbitrary boundary conditions [78,80,90], once the corresponding DS formulation has been developed. The procedure is simply replacing the DS formulations of Equation (4) by the formulations in [78,80,89,90]. Likewise, the rod/beam formulations could be quite general to consider composite ones as well. Moreover, this type of point connections have very important practical applications as mentioned earlier in the Introduction, and whose DS formulation can be accomplished in a concise manner.

The Global DS Formulation and Mode Count for Line Nodes Subjected to Point Supports/Connections
We provide the formulation of global DS and extension of the W-W algorithm for rigid point supports/connections (Section 3.2.1), elastic point supports/connections without slave degrees of freedom (Section 3.2.2) or connected to substructures with slave degrees of freedom (Section 3.2.3).

Rigid Point Supports/Connections
When a plate is subjected to rigid point supports/connections as shown in Figure 6.
Putting Equations (33) and (35) together leads to the following matrix relationship: There are two different situations in terms of Lagrange multipliers: (i) rigid point supports; and (ii) rigid point coupling connection. The only difference is in the expression of E P . The usage of Lagrange multiplier in the DSM will introduce extra DOFs, in this case, the sign count of these extra DOFs, which should be eliminated from the system under consideration. Assuming that there are N P rigid point constraints, the associated mode count formula will be as follows: where, J 0plm is the number of eigenvalues for the plate element, which would still lie between 0 and ω * and N P is the number of rigid point supports.
In Equation (37), we have Therefore, Equation (37) could also take the form which is the previously extended W-W algorithm proposed by Williams and Anderson [38,84]. It worth emphasizing that special care should be taken when applying rigid point coupling constraints to avoid overdetermined constraints; otherwise, the above mode count algorithm will lead to misleading results.

Elastic Point Supports/Connections without Slave Degrees of Freedom
When the parent structure is connected to elastic point supports/connections without slave degrees of freedom, such as some representative cases shown in Figure 7, the equilibrium conditions at these point supports/connections can be written as where K P is the DS matrix of the elastic supports or connections, such as spring supports, rods and beams with the other end clamped. By combining Equations (33), (35) and (40) together, we can obtain the DS formulation of line nodes connected to point elastic supports as which is applicable to a wide range of parent structures with line nodes connected to point nodes without slave degrees of freedom, see Figure 7. Moreover, the corresponding mode count formula takes the following unique form: where, J 0plm and J 0m p are the number of eigenvalues for plate and beam/bar element, which would still lie between 0 and ω * , respectively. N P is the degree of point connection freedom. Note that E P takes different forms for each case. Notably, K P in Equation (41) can represent three cases.

1.
Parent structure subjected to massless spring point supports (Figure 7a). Assuming that the DS of the point supports are k p . In Equation (40), the matrix formulation K P can be written as K P = diag[k p ], (p ∈ [1, P]).

2.
Parent structure subjected to a rod or beam point connection with the other end clamped (Figure 7b). The formulation for K P will then be the same as Equation (40) but K P represents the supported/connected DS matrix of beam or rod.

3.
Parent structure with point elastic coupling constraints (Figure 7c). In this case, if two plates are coupled by springs, K P can be described by Equation (20). However, when two plates are coupled by a rod or beam, K P represents the connected DS matrix of beam or rod.

Elastic Point Supports/Point Connections to Substructures with Slave Degrees of Freedom
The substructure is connected with the parent structure (DS matrix K P ) with the d P whilst other DOFs are related to displacement d s . A parent structure can be connected with any number of substructures through point connections, where the substructures have point nodes either in the form of 1D analytical DS elements or 1D or 2D numerical finite element elements, see Figure 8. The DS formulation of the substructures can be rewritten in the form f sub = K sub d sub (43) i.e., f p and d p are generalized forces and displacements at point node/connections to plate, f s and d s denote the generalized forces and displacements of slave nodes without connecting to plates, K sub is the DS matrix of the substructure. We can then formulate the global DS matrix of point connections to other substructures by putting Equations (33), (35) and (44) together.
Moreover, the corresponding mode count formula take the following unique form where, J 0plm and J 0m sub are the number of eigenvalues for plate and beam/bar element which would still lie between 0 and ω * , respectively. N P is the DOF of the connection point between the plate and the substructure. It is easily seem that the case formulated by Equation (45) is a general case of those formulated by Equation (41), which can be implemented to the computer program covering all possible cases.
The solution procedure of the proposed dynamic stiffness formulation as well as the associated W-W algorithm for the free vibration analysis of built-up structures is demonstrated by the flow chart given in Figure 9.

Results and Applications
The new DS formulations and the corresponding enhanced Wittrick-Williams algorithm have been implemented in a Matlab code. The DS models of several representative built-up structures are constructed and computed results are compared and contrasted with those computed by commercial FE software ANSYS. For representative boundary conditions, a simply supported edge is denoted by the letter S and a free edge is represented by the letter F. In the examples of this paper, unless otherwise specified, the size of the Levy-type plate is a = b = 1 m, the boundary condition for the typical FSFS is shown in Figure 10, Young's modulus is taken as E pl = 7.2 × 10 10 Pa, Poisson ratio ν pl = 0.3 and the mesh size of plate FEM elements is chosen to be 5 mm.

Example 1: Plate with Rigid Point Supports
In the numerical implementation of the DSM theory, different terms of Levy solutions are superposed to model a plate element. It is important to examine the convergence rate of the method and to verify the DS formulation of Equation (36) and the extension of W-W algorithm of Equation (37). Using the classical thin plate theory, we study the convergence and computational efficiency of a plate supported by three rigid points at (0,0.25,0), (0.4,0.5,0), (1,0.75,0) as shown in Figure 11. The plate parameters of density ρ pl = 2700 kg/m 3 , thickness h pl = 0.002 m and other parameters have been given previously. The FE mesh is chosen to be 100 × 100, 200 × 200 and 500 × 500, respectively. It is clearly demonstrated from Table 1 that the current DSM converges very rapidly to exact solutions. In addition, the first ten natural frequencies obtained by the DSM match very well with the FE results, the difference might be due to the discrepancy of classical plate theory adopted in the DSM and the advanced plate theory adopted in the FEM, and the calculation time is obviously better than commercial software ANSYS. When the half-wavenumber m = 25, the proposed DSM is 23 times faster than ANSYS with a 200 × 200 mesh. It is shown that the DSM for plate supported by rigid points has excellent convergence rate and computational efficiency. Considering the accuracy of the calculation performance and efficiency, the following calculation examples take m = 25.

Example 2: Plates Supported at Points by Massless Springs
In this example, we studied a plate subjected to three elastic point supports to verify the DS formulation of Equation (41) and the extension of the W-W algorithm of Equation (42). As illustrated in Figure 12, the coordinates of three spring points supporting the plate are (0,0.25,0), (0.4,0.5,0) and (1,0.75,0), the plate is formed by two plate elements and other parameters of the plate are the same as Example 1. Table 1. The first 10 natural frequencies for a plate supported by three rigid points as shown in Figure 11, by the current DSM with different m, compared with FEM solutions with three different meshes.  We compute the first three natural frequencies of the structure with a range spring stiffness covering K s = 0∼2 × 10 3 kN/m as shown in Figure 13. It can be seen that the first three natural frequencies of the structure approach the results of plate subject to rigid point supports as shown in Table 1. Next, the 1st∼10th , 20th, 30th, 40th and 50th natural frequencies under the three massless springs with different stiffnesses of K s = 10 N/m, 3 kN/m and 20 kN/m are tabulated in Table 2. The 1st, 2nd, 8th and 9th natural mode shapes with K s = 3 kN/m computed by both DSM and FE package ANSYS are compared and contrasted in Figure 14. Table 2. The first 10 natural frequencies and the 20th, 30th, 40th and 50th natural frequencies for a plate subjected to three elastic spring point supports as shown in Figure 12.   As can be seen in Table 2, the DSM results are in good agreement with the FE results. When the spring stiffness K s increases, the first five natural frequencies increase as expected, but the effect on the natural frequencies after the 10th mode is much less significant. This is not surprising since the boundary conditions have a much smaller effect on the vibration behaviors in the high frequency range than in the low frequency range because the characteristics of overall modes are very different from local modes. In addition, four representative mode shapes of the structures computed by both the DSM and the FEM are presented and contrasted in Figure 14, which demonstrates good agreement.

Example 3: A Plate Connected by a Rod and/or a Beam
In this example, we examine the free vibration behaviors of a plate connected to a rod and/or a beam, respectively, to verify the DS formulation of Equation (45) and the extension of the W-W algorithm of Equation (46). As is shown in Figure 15, a rod is perpendicularly connected to the parent plate and the other end of the rod is free see Figure 15a, and a beam horizontally connected to the plate and the other end of the beam is rotated by θ, i.e., the beam is constrained see Figure 15b. Both the rod and the beam are connected to the plate at position (0,0.25,0). We computed from the 1st to 10th natural frequency for the plate connected to the rod and/or the beam, as shown in the left and right half of Table 3, respectively. It can be concluded from the result of Table 3 that the first 10 natural frequencies computed by the DSM and FE match very well with each other for two models described in Figure 15. Four representative mode shapes of the model described in Figure 15b are shown in Figure 16. Comparisons are made for both DSM and FE results, which show excellent agreements. Table 3. The first 10 natural frequencies for the plate connected to a rod and/or a beam shown in Figure 15.

Example 4: Two Plates Coupled by a Spring or a Beam
This example can be applied to verify the DS formulation of Equation (41) and the extended W-W algorithm of Equation (42). As shown in Figure 17a, two plates coupled by a spring or a beam: the spring is perpendiculariy connected to two parallel plates at position (0,0.25,0.3) and (0,0.25,0). The plate thickness h pl = 0.005 m with density ρ pl = 2700 kg/m 3 and the other parameters are the same as Example 1, the spring stiffness K s = 1.2 kN/m. As shown in Figure 17b for the next example, a beam is horizontally connected to two plates at (1, 0.25, 0) and (1.4, 0.25, 0). The plate thickness h pl = 0.006 m, density ρ pl = 2700 kg/m 3 ; the beam Young's modulus E r = 2 × 10 11 Pa density ρ r = 7800 kg/m 3 , cross section B b = H b = 0.004 m, length l b = 0.4 m and rotation angle θ constrained. Next, we apply both the proposed DSM and the FE package ANSYS to compute the 1st∼10th natural frequencies (there are shown in Table 4) and four representative mode shapes (shown in Figures 18 and 19) for the two cases shown in Figure 17, and contrast the four representative mode shapes computed by both DSM and FE package ANSYS. Table 4. The first 10 natural frequencies for two plates coupled by a spring with K s = 1.2 kN/m shown in Figure 17a and by a beam shown in Figure 17b. It can be ascertained from Table 4 that the first ten natural frequencies calculated by the DSM and FEM match extremely well with each other for the two cases described in Figure 17. Four representative mode shapes of the two cases computed by the DSM and the FEM are compared in Figures 18 and 19, which show good agreements.

Conclusions
The DSM and the Wittrick-Williams algorithm are enhanced and subsequently used to model 2D line nodes with 1D point node supports/connections built-up structure accurately and efficiently. Motivated by the work of Williams and Anderson [84], the proposed method is a unified approach towards universal DS modeling and a further extension for built-up structures with line nodes connected to point nodes. The proposed developments have the following advantages: 1.
The theory is systematic and general with a wider scope for applications. The W-W algorithm is extended to cover a wide range of cases, including plates with point rigid constraint, point elastic support, rod support, rigid support, beam constraint, point elastic coupling constraint, plate-rod coupling, and plate-beam coupling.

2.
The modeling procedure is more direct and convenient with much better numerical stability. The work by Williams and Anderson [84] are based on dynamic flexibility formulation, which involved matrix inversions that might become numerically unstable for larger systems with a larger number of DOFs.
The effect of the 2D line nodes supported/connected by 1D point nodes on the natural frequencies and mode shapes of the built-up structure are demonstrated by numerical results which show excellent agreement with ANSYS results. The proposed method is computationally efficient and gives highly accurate results, which can be further extended to complex structural systems.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Acknowledgments:
We would like to thank the reviewers for their valuable comments and suggestions, which helped us to improve the quality of the article.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. DS Matrix for Classical Plate Theory
The 4 × 4 DS matrix of Equation (7) including six independent terms k 11 , k 12 , k 13 , k 14 , k 22 , k 24 . Thus, K pl (ω) can be expressed as: h2 − R 2 S h1 S h2 k 13 = (L 2 − L 1 )(r 2 C h2 S h1 − r 1 C h1 S h2 ) k 14 = (r 2 R 1 − r 1 R 2 )(r 2 S h2 − r 1 S h1 ) k 22 = (C h2 − C h1 )(r 2 R 1 − r 1 R 2 ) k 24 = (L 1 − L 2 )(r 2 S h1 − r 1 S h2 ) For the Wittrick-Williams algorithm, the J 0 of the plate element problem is solved by the indirect method [77]. According to the Wittrick-Williams algorithm, the mode count J sm of the plate element with all nodal boundaries simply supported when the halfwavenumber in the y direction is m can be given by equation Equation (21), which can be recast as where J 0m is J 0 when the half-wavenumber in the y direction is m and K s plm (ω * ) is the DS matrix K plm (ω * ) for a plate element with all nodal boundaries simply supported for a given half-wavenumber m.

Appendix B. DS Matrix for the Classical Rod Thory
The dynamic stiffness formula for axial vibration of classical rod can be obtained, namely Equation (12): where e 1 = EA l µ cot µ, e 2 = − EA l µ csc µ, µ = ρAω 2 l 2 EA The W-W algorithm solving the natural frequencies of the rod element in the same form as equation Equation (21), and J 0r of the rod element expressed as follows: where I NT is the integer function.