Extended analytical strip method for free vibration analysis of rectangular thin plates with arbitrarily varying parameters along one edge

Abstract

Thin plates with properties varying along a single in-plane axis arise in many engineering applications, yet their analysis is often hindered by variable-coefficient governing equations that challenge existing analytical methods. This study presents an extension of the analytical strip method (ASM) for the free vibration analysis of rectangular thin plates with varying material properties and thickness along one edge. Unlike previous ASM formulations that treat only stepped plates and homogeneous laminates, the proposed approach can model arbitrarily varying thickness, density, and material properties. The governing differential equation is derived under the local Kirchhoff–Love assumptions, reduced via separation of variables and employing Lévy-type solution in the simply supported direction, and is discretized into strips along the varying axis. Within each strip, exponential trial functions yield closed-form homogeneous solutions. Continuity conditions across inter-strip boundaries enforce kinematic and load smoothness; while clamped, simply supported, or free boundary conditions (or combinations thereof) apply at the two outer longitudinal edges of the boundary strips. Three numerical examples are presented in order to validate the method. One example is chosen from the literature and the others are designed so as to have analytical solutions despite strong nonlinearities. The results, benchmarked against finite element and analytical solutions, demonstrate rapid convergence and excellent accuracy for plates with variable material properties and thickness, even with relatively much smaller number of elements.

1. Introduction

Global advances in lightweight structures, aerospace components, wind-turbine skins, and civil-engineering panels increasingly rely on thin plates, stiffness, thickness, and density of which vary intentionally along one direction to improve efficiency, reduce weight, or tune vibrational performance (Hassan and Saeed, 2024; Madsen et al., 2012; Sonmez, 2017; Stanford and Jutte, 2014; Wang et al., 2019). This trend is also reflected in the vibration literature, where one-direction thickness tapering and in-plane material non-homogeneity have been shown to alter natural frequencies and mode shapes of thin plates (Lal and Saini, 2017; Manna, 2012). Because even small changes in these properties may shift resonant frequencies and alter dynamic response, accurately predicting the vibration behavior of such non-uniform plates remains essential for reliable structural design.

The analytical strip method (ASM), originally formulated by Harik and Pashanasangi (1985) and later extended to thin rectangular orthotropic plates by Harik and Salamoun (1986), then to antisymmetric laminates by Sun (2009) and Sun and Harik (2010), provides closed-form single-series solutions for rectangular plates with two opposite simply supported edges subjected to patch, line, and concentrated loads. However, these works assume homogeneous or thickness-stepped plates/laminates with constant cross-sectional properties along one axis.

Over the past few decades, ASM has been successfully extended and applied to a variety of problems in plate and shell analysis. Harik and Salamoun (1988) modeled stiffened plates analytically by idealizing them as strips connected to beam segments (utilized as stiffeners). Harik and Andrade (1989) applied the method to elastic stability analysis of plates with stepped thickness. Harik et al. (1992) addressed free vibration of orthotropic plates using ASM. Sun and Harik (2015) extended the formulation to stiffened and continuous antisymmetric laminates. ASM has also been applied to shell problems (Perkins, 2017; Perkins and Harik, 2017) and to curved plates and bridge decks (Harik and Pashanasangi, 1985). These developments confirm the flexibility and accuracy of ASM in handling patch, uniform, and point loads, as well as varying support conditions and stiffeners.

Modern engineering applications increasingly employ functionally graded and geometrically tapered plates (e.g., wing skins, bridge decks, ship hulls) with continuously varying properties to optimize performance and weight (Harik and Pashanasangi, 1985). Numerical methods such as finite elements handle such non-uniformities but at the cost of mesh generation and potential numerical dispersion. An analytical strip framework for continuously varying plates would combine high accuracy with conceptual simplicity (Harik and Andrade, 1989; Harik et al., 1992; Harik and Salamoun, 1986, 1988).

For such non-uniform structures, finite element methods remain the dominant general-purpose tool, with recent p-version and higher-order formulations reported for porous and functionally graded plates, shells, and beam-type members, including porous and carbon nanotube (CNT) reinforced sandwich plates, bi-directional porous shells, thermally loaded porous functionally graded (FG) material beams, porous nanocomposite beams on viscoelastic foundations, variable-cross-section sandwich beams, beams resting on Winkler–Pasternak foundations, profile optimization of beams with tip mass for energy harvesting, and axially FG curved beams on orthotropic Pasternak foundation (Bentrar et al., 2023; Bousmaha et al., 2026; Lakhdar et al., 2024; Youzera et al., 2025; Tounsi et al., 2026; Belabed et al., 2024e, 2024a, 2024b, 2024c, 2024d; Karadag et al., 2021; Ermis et al., 2022). Related finite element developments in structural mechanics more broadly include homogenization-based nonlinear formulations for masonry walls (Meftah et al., 2024).

At the same time, semi-analytical assembly-type methods continue to be attractive in one-dimensional structural vibration because they preserve exact local homogeneous solutions while assembling only boundary and interface conditions numerically. A related example is the numerical assembly technique/method (NAT/NAM), developed for beam and frame vibration problems, where high accuracy is obtained without conventional spatial interpolation of the field variable (Kramer and Gfrerer, 2024). More broadly, recent vibration studies have also pursued analytical and semi-analytical treatments of non-uniform members and resonant systems, including post-critical nonlocal strain-gradient beams, hemispherical resonators with mass defects, and acoustic-black-hole beam configurations (Alam et al., 2025; Gao et al., 2021; Liu et al., 2025).

Compared with the finite element method (FEM), the present strip formulation does not approximate the field equations within each strip by low-order interpolation functions; instead, it satisfies the local strip equation exactly after coefficient freezing. Compared with Ritz, Galerkin, or spectral approaches, the present method enforces boundary and inter-strip conditions directly in physical space, which makes it convenient for mixed edge conditions and for extracting multiple modes once the frequency search interval is prescribed. In this sense, the extended ASM should be viewed as a semi-analytical bridge between exact Lévy-type plate solutions and fully discretized numerical models.

This study extends the ASM to plates exhibiting variations in thickness, material properties such as elastic moduli and Poisson’s ratios, as well as mass density along one axis (y). The extension retains the classical separation of variables in the other direction (x, simply supported direction) while modeling y-variations by subdividing the plate into narrow strips and solving fourth-order ODEs with constant coefficients within each strip using exponential functions. Appropriate matching of displacements and their derivatives up to the third order ensures necessary and required smoothness across strip boundaries. The result is a semi-analytic solution valid for the free vibration of rectangular thin plates where the axis of variation (in material properties like modulus or density) is orthogonal to a simply supported direction. The classical ASM is generalized from homogeneous and stepped-property plates to plates with continuously as well as discretely varying thickness, density, and in-plane stiffness parameters in one dimension. This is achieved without abandoning the semi-analytical character of the method: after the Lévy-type reduction, each strip is still solved in closed form, and only the strip coupling and eigenvalue extraction remain numerical. Thus, the method fills the gap between exact Lévy-type plate solutions, which require more restrictive coefficient structures, and fully discretized finite element models.

For verification, three examples are presented, two of which allow analytical solutions despite nonlinearly varying properties and thickness. One of the examples was taken from the literature; Kang and Kim (2008) and Singh and Saxena (1996). In all cases, ASM results were compared to the results from the literature or those obtained via analytical solutions or finite element analysis (ANSYS), whichever were applicable.

It was observed that ASM was superior to other methods and quickly converged to analytical solutions in predicting natural frequencies and mode shapes of thin plates with variable properties and thickness.

2. Method

Moments (see Figure 1) integrated across the thickness of an infinitesimal plate element in terms of the stresses can be written as (Graff, 1975)

$\left[\begin{array}{c}\hfill M_x\hfill \\ \hfill M_y\hfill \\ \hfill M_{xy}\hfill \end{array}\right]={\int }_{-\frac{h}{2}}^{\frac{h}{2}}\left[\begin{array}{c}\hfill {\sigma }_x\hfill \\ \hfill {\sigma }_y\hfill \\ \hfill {\sigma }_{xy}\hfill \end{array}\right]zdz$

(1)

Figure 1.

Infinitesimal thin plate element with moment loads. Normal and shear forces are not shown.

where z is the transverse direction (with respect to the mid-surface) and h is the thickness of the thin plate and the stresses in relation to the strains for the plane stress case are (Reddy, 2006)

$\left[\begin{array}{c}\hfill {\sigma }_x\hfill \\ \hfill {\sigma }_y\hfill \\ \hfill {\sigma }_{xy}\hfill \end{array}\right]=\left[\begin{array}{ccc}\hfill Q_{11}\hfill & \hfill Q_{12}\hfill & \hfill 0\hfill \\ \hfill Q_{12}\hfill & \hfill Q_{22}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill Q_{66}\hfill \end{array}\right]\left[\begin{array}{c}\hfill {ϵ}_x\hfill \\ \hfill {ϵ}_y\hfill \\ \hfill {ϵ}_{xy}\hfill \end{array}\right]$

(2)

where Q _rs are the plane stress-reduced stiffnesses and are obtained from:

$\left[\begin{array}{c}\hfill {ϵ}_x\hfill \\ \hfill {ϵ}_y\hfill \\ \hfill {ϵ}_z\hfill \\ \hfill {ϵ}_{yz}\hfill \\ \hfill {ϵ}_{zx}\hfill \\ \hfill {ϵ}_{xy}\hfill \end{array}\right]=\left[\begin{array}{cccccc}\hfill S_{11}\hfill & \hfill S_{12}\hfill & \hfill S_{13}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill S_{12}\hfill & \hfill S_{22}\hfill & \hfill S_{23}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill S_{13}\hfill & \hfill S_{23}\hfill & \hfill S_{33}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill S_{44}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill S_{55}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill S_{66}\hfill \end{array}\right]\left[\begin{array}{c}\hfill {\sigma }_x\hfill \\ \hfill {\sigma }_y\hfill \\ \hfill {\sigma }_z\hfill \\ \hfill {\sigma }_{yz}\hfill \\ \hfill {\sigma }_{zx}\hfill \\ \hfill {\sigma }_{xy}\hfill \end{array}\right]$

(3)

where for the orthotropic thin plate model (plane-stress case) σ _z = ϵ _yz = ϵ _zx = 0, in-plane shear is neglected σ _yz ≃ σ _zx ≃ 0.

Inverse of the reduced compliance matrix in equation (3) yields the plane stress-reduced stiffnesses Q _rs . For the isotropic and homogeneous case, the stiffnesses become Q₁₁ = Q₂₂ = Q₁₂/ν = Q₆₆/(1 − ν) = E/(1 − ν²) and h (where E is Young’s Modulus and ν is Poisson’s ratio) constant, the reduced stiffness matrix elements become those of the isotropic and homogeneous case. Strains are functions of displacements, which in this case involve only transverse displacement w.

Combining equation (1) and the non-trivial (one force and two moment) equations of motion along with stress-strain relations in equation (2), a single equation (equation (4)) in terms of w is obtained. As a result, the transverse displacement w(x, y) of a thin plate with spatially variable coefficients along y-axis (a _i = a _i (y), i = 1, 2, …, 8) satisfies

$\begin{array}{ccc}\hfill a_1\frac{{\partial }^{4}w}{\partial x^4}+a_2\frac{{\partial }^{2}w}{\partial x^2}+a_3\frac{{\partial }^{3}w}{\partial x^2\partial y}+a_4\frac{{\partial }^{4}w}{\partial x^2\partial y^2}& \hfill & \hfill \\ \hfill +a_5\frac{{\partial }^{2}w}{\partial y^2}+a_6\frac{{\partial }^{3}w}{\partial y^3}+a_7\frac{{\partial }^{4}w}{\partial y^4}+q+a_8\frac{{\partial }^{2}w}{\partial t^2}=0& \hfill & \hfill \end{array}$

(4)

where q = q(x, y) is the applied load distribution. The coefficients a _i (y) are given as follows:

$\begin{array}{ll}\hfill a_1& =Q_{11}\hfill \end{array}$

(5a)

$\begin{array}{ll}\hfill a_2& =\left(\frac{6{(h^{\prime })}^2}{h^2}+\frac{3h^{\prime \prime }}{h}\right)Q_{12}+\frac{6h^{\prime }}{h}{Q}_{12}^{\prime }+{Q^{\prime \prime }}_{12}\hfill \end{array}$

(5b)

$\begin{array}{ll}\hfill a_3& =\frac{6h^{\prime }}{h}{Q}_{12}+2Q_{12}^{\prime }+\frac{6h^{\prime }}{h}{Q}_{66}+2Q_{66}^{\prime }\hfill \end{array}$

(5c)

$\begin{array}{ll}\hfill a_4& =2\left(Q_{12}+Q_{66}\right)\hfill \end{array}$

(5d)

$\begin{array}{ll}\hfill a_5& =\left(\frac{6{(h^{\prime })}^2}{h^2}+\frac{3h^{\prime \prime }}{h}\right)Q_{22}+\frac{6h^{\prime }}{h}{Q}_{22}^{\prime }+{Q^{\prime \prime }}_{22}\hfill \end{array}$

(5e)

$\begin{array}{ll}\hfill a_6& =\frac{6h^{\prime }}{h}{Q}_{22}+2Q_{22}^{\prime }\hfill \end{array}$

(5f)

$\begin{array}{ll}\hfill a_7& =Q_{22}\hfill \end{array}$

(5g)

$\begin{array}{ll}\hfill a_8& =12\frac{\rho }{h^2}\hfill \end{array}$

(5h)

where ρ = ρ(y) is the density of the plate, and (…)′ denotes the derivative with respect to y.

Assuming simply supported edges at x-boundaries and strip discretization (see Figure 2) w can be given in series form as follows:

$w(x,y,t)=\sum _{l=1}^{\infty }\sum _{k=1}^{\infty }{Y}_{kl}(y)\sin (k\pi x/L_x)\sin ({\omega }_{kl}t+\varphi )$

(6)

where L _x is the x-axis length of the strips and k and l are the x and y mode indices, respectively. Substituting into equation (4) and exploiting orthogonality in x reduces the problem to a family of fourth-order ODEs in y for each Y _kl , simplified as Y:

$\begin{array}{ccc}\hfill a_7{Y}^{\prime \prime \prime \prime }+a_6{Y}^{\prime \prime \prime }+(a_5-a_4{\mathrm{\Lambda }}_x)Y^{\prime \prime }-a_3{\mathrm{\Lambda }}_x{Y}^{\prime }& \hfill & \hfill \\ \hfill +(a_1{\mathrm{\Lambda }}_x^2-a_2{\mathrm{\Lambda }}_x-a_8{\omega }^2)Y=0& \hfill & \hfill \end{array}$

(7)

where ${\mathrm{\Lambda }}_x=\frac{k^2{\pi }^2}{L_x^2}$ . The domain where the solution of the ODE is searched is discretized into N − 1 strips. The functions a _i (y) are only assumed to be finite and integrable in each strip.

Figure 2.

Segmented plate with arbitrary thickness along y-axis.

The values of a _i (y) for each segment is calculated as the mean of the functions $a_{ni}=\frac{1}{y_{n+1}-y_n}{\int }_{y_n}^{y_{n+1}}{a}_i(y)dy$ , where y _n and y_n+1 are the y-axis boundaries of the n ^th segment. In general, any continuous function or functions with finitely many finite discontinuities are admissible. If the functions a _i (y) are not changing rapidly and are continuous, values of the functions at midpoints of segments could also be expected to yield approximate solutions. In this way, for each segment the corresponding constant coefficient ODE becomes

$c_{n4}{Y}_n^{\prime \prime \prime \prime }+c_{n3}{Y}_n^{\prime \prime \prime }+c_{n2}{Y}_n^{\prime \prime }+c_{n1}{Y}_n^{\prime }+c_{n0}{Y}_n=0$

(8)

where the coefficients c _nj incorporate a _i = a _i (y), i = 1, 2, …, 8 and j = 0, 1, …, 4.

Since each equation is a constant coefficient fourth-order linear ODE, solutions are

$Y_n=A_{n1}{e}^{{\lambda }_{n1}y}+A_{n2}{e}^{{\lambda }_{n2}y}+A_{n3}{e}^{{\lambda }_{n3}y}+A_{n4}{e}^{{\lambda }_{n4}y}$

(9)

$n=\mathrm{1,2},\dots ,N-1$

(10)

where λ _nj ’s are the zeros of the characteristic equation (roots of the characteristic polynomial) that can be real or complex and are functions of ω.

At strip interfaces, the following continuity conditions, providing sufficient smoothness, are enforced in case of continuous parameter functions:

$\begin{array}{l}\hfill Y_n(y_{n+1})=Y_{n+1}(y_{n+1})\end{array}$

(11a)

$\begin{array}{l}\hfill Y_n^{\prime }(y_{n+1})=Y_{n+1}^{\prime }(y_{n+1})\end{array}$

(11b)

$\begin{array}{l}\hfill {Y^{\prime \prime }}_n(y_{n+1})={Y^{\prime \prime }}_{n+1}(y_{n+1})\end{array}$

(11c)

$\begin{array}{l}\hfill Y_n^{\prime \prime \prime }(y_{n+1})=Y_{n+1}^{\prime \prime \prime }(y_{n+1})\end{array}$

(11d)

If there are discontinuities in the parameter functions, then instead of third and fourth-order derivative continuity, the internal shear and bending moment continuity conditions should be used since there are no external loads. The shear and bending moment functions can be given as:

$\begin{array}{ccc}\hfill M_y=\frac{h^3}{12}{K}_{yy}X{Y}^{\prime \prime }-{\mathrm{\Lambda }}_X\frac{h^3}{12}{K}_{xy}XY& \hfill & \hfill \end{array}$

(12)

$\begin{array}{ccc}\hfill V_y=-\frac{h^2}{12}X\left[h{K}_{yy}{Y}^{\prime \prime \prime }+(h{K}_{yy}^{\prime }+3h^{\prime }{K}_{yy})Y^{\prime \prime }\right.& \hfill & \hfill \\ \hfill \left.-(h{K}_{xy}+2h{D}_{xy}){\mathrm{\Lambda }}_x{Y}^{\prime }-(h{K}_{xy}^{\prime }+3h^{\prime }{K}_{xy}){\mathrm{\Lambda }}_{x}Y\right]& \hfill & \hfill \end{array}$

(13)

In general, a total of 4N equations are obtained from the y-axis boundary and continuity conditions. These equations can be given in matrix form involving a 4N × 4N coefficient matrix (representing the boundary and continuity conditions) $\tilde{M}$ and a 4N × 1 coefficient vector $\overline{A}$ , such that $\tilde{M}\overline{A}=\overline{0}$ , where $\tilde{M}$ is a function of ω. For non-trivial solutions, one must have $\det (\tilde{M}(\omega ))$ = 0.

An example 8 × 8 matrix formed from a two element (N = 3) ASM application is shown in Figure 3, where all of the λ’s are functions of ω. In general, $\tilde{M}$ is a block diagonal matrix.

Figure 3.

Coefficient matrix ${\tilde{M}}_{8\times 8}$ of a two element (N = 3) ASM application example with two edges clamped (SCSC), where λ′s are functions of ω.

Due to the fact that $\tilde{M}$ is usually a large matrix, elements of which are nonlinear functions of ω, solving the determinant equation by the usual methods becomes practically unfeasible. Instead, in a preliminary frequency sweep, the first vanishing eigenvalue is determined roughly. This amounts to solving the following equation:

$min(abs(eig(\tilde{M}(\omega ))))=0$

(14)

A vanishing eigenvalue is a sufficient condition for having a vanishing determinant. Equation (14) is evaluated over a prescribed frequency interval by a preliminary sweep between 10⁻⁴ and 1000 Hz. Some frequencies are identified as local minima candidates. If the algorithm does not find a candidate in any run, the resolution is first increased, then the search domain is increased until a candidate is found. Each candidate is then refined by a local one-dimensional minimization procedure, after which the corresponding mode shape is reconstructed and checked against the boundary residuals. In order to increase the accuracy, the result is improved by minimizing the maximum of the absolute values of the normalized deflections Y(y) and/or their relevant derivatives, which are required to vanish at the boundaries. See Figure 4 for the combinations of boundary conditions (BC) used in the numerical studies. With each result, the search domain is updated accordingly. If the algorithm finds nonconsecutive modes in a search interval, then attempts to find the missing modes by shortening the interval to that between these nonconsecutive modes.

Figure 4.

Combinations of boundary conditions present in the analyses.

In case of eigenvalue multiplicity, the algorithm determines eigenmodes correctly if the first (x) mode index numbers are different. When the first mode index is the same, the current algorithm is not designed to resolve the corresponding mode shapes due to the decreased rank 4N − s of $\tilde{M}$ , instead of 4N − 1, where s > 1 is the degree of multiplicity. However, the algorithm can be expanded in a straightforward fashion to investigate such cases, albeit being rare.

The present formulation is intended for thin plates for which Kirchhoff–Love assumptions remain appropriate; that is, transverse shear deformation and rotary inertia are assumed negligible compared to bending effects. The method is therefore most suitable when the local thickness remains small relative to the in-plane dimensions, even if that thickness varies appreciably along one direction. For moderately thick plates or cases involving strong local shear effects, a Mindlin–Reissner extension would be more appropriate. In addition, because a Lévy-type separation is employed, the present development is restricted to plate configurations having two opposite simply supported edges.

A flowchart outlining the solution procedure is presented in Appendix A.

3. Results

In this section, three examples are presented in order to demonstrate the validity and performance of the extended ASM method developed in this study. The first example is taken from a previous study (Kang and Kim, 2008), the remaining two are designed so as to serve as benchmarks for the method.

3.1. Example 1

Consider a thin rectangular plate (material: E = 200 GPa, ν = 0.3, ρ = 7850 kg/m³) with a linearly varying thickness with all edges simply-supported. Two cases of plate thickness, one following a linear profile, h(y) = 0.001(1 + 0.6y), the other following a quadratic profile h(y) = 0.001(1 + 0.6y²) are studied in Kang and Kim (2008). Both gradually thin toward one end in the domain x ∈ [0, 0.5] and y ∈ [0, 1].

Table 1 presents the results obtained by MATLAB for ASM, ANSYS, and transfer matrix method (TMM) of Kang and Kim (2008) and Rayleigh-Ritz (RR) of Singh and Saxena (1996). The ASM solution introduced in this study, the extended ASM, with even relatively few strips (e.g., 80 strips) shows excellent agreement with a fine-mesh finite element analysis (17956 elements in ANSYS) and with published results (Kang and Kim, 2008), demonstrating the method’s accuracy and rapid convergence (Table 2) for non-uniform thickness.

Table 1.

Natural frequencies (Hz) of an SSSS rectangular plate with thickness varying along y: linear and quadratic tapers. Comparison of ASM (80 strips), TMM (Kang and Kim, 2008; 40 terms for linear, 30 for quadratic), RR (Singh and Saxena, 1996), and ANSYS (17956 elements).

Thickness variation	Method/Ref.	Mode shape
		(1,1)	(1,2)	(1,3)	(2,1)	(2,2)	(2,3)
Linear	ANSYS	15.36	24.74	40.11	50.80	62.15	77.47
	TMM	15.40	24.79	40.16	50.91	62.33	∗
	RR	15.37	24.77	40.28	∗	∗	∗
	ASM	15.37	24.76	40.14	50.81	62.18	77.53
Quadratic	ANSYS	14.04	22.82	36.96	45.96	56.91	71.40
	TMM	13.84	22.79	37.13	45.94	57.07	∗
	ASM	14.05	22.83	36.98	45.97	56.94	71.45

Table 2.

Convergence for the linearly tapered SSSS plate (Kang and Kim, 2008): ASM frequencies (Hz) versus number of strips; ANSYS frequencies (Hz) versus mesh refinement.

BC	Method	Elem. #	Mode shape
			(1,1)	(1,2)	(1,3)	(2,1)	(2,2)	(2,3)
SSSS	ASM	2	15.08	25.01	39.62	49.72	63.27	76.23
		3	15.26	24.53	40.52	50.30	61.61	78.66
		4	15.31	24.65	39.88	50.56	61.86	77.06
		5	15.33	24.70	40.01	50.67	62.01	77.25
		6	15.34	24.72	40.06	50.72	62.07	77.37
		7	15.35	24.73	40.09	50.74	62.11	77.43
		8	15.35	24.74	40.10	50.76	62.13	77.46
		9	15.36	24.75	40.11	50.77	62.14	77.48
		10	15.36	24.75	40.11	50.78	62.15	77.49
		11	15.36	24.75	40.12	50.79	62.15	77.50
	ANSYS	144	15.34	24.78	40.49	51.32	63.96	83.82
		182	15.34	24.75	40.26	51.04	62.99	80.50
		240	15.35	24.73	40.14	50.91	62.50	78.67
		676	15.36	24.74	40.10	50.81	62.17	77.52
		2550	15.36	24.74	40.11	50.80	62.15	77.48
		10100	15.36	24.74	40.11	50.80	62.15	77.47
		17956	15.36	24.74	40.11	50.80	62.15	77.47

3.2. Example 2

In this example, a square plate, with two opposite edges simply supported and the other two remaining edges being either clamped (C), simply-supported (S), or free (F), is considered. The thickness and density are allowed to vary exponentially along the y direction as $h(y)=c_h{e}^{{\lambda }_{h}y}$ and $\rho (y)=c_{\rho }{e}^{{\lambda }_{\rho }y}$ , respectively, with the parameters related as λ _ρ = 2λ _h . This particular choice of parameters causes the y-dependent terms in the plate’s ODE (after separation of variables) to simplify considerably. Further, if the material is taken as isotropic with constant elastic moduli yields the following fourth-order ordinary differential equation in y with constant coefficients.

$\begin{array}{ccc}\hfill Q_{22}{Y}^{\prime \prime \prime \prime }+6{\lambda }_h{Q}_{22}{Y}^{\prime \prime \prime }& \hfill & \hfill \\ \hfill +(9{\lambda }_h^2{Q}_{22}-2{\mathrm{\Lambda }}_x(Q_{12}+Q_{66}))Y^{\prime \prime }& \hfill & \hfill \\ \hfill -6{\lambda }_h{\mathrm{\Lambda }}_x(Q_{12}+Q_{66})Y^{\prime }& \hfill & \hfill \\ \hfill +({\mathrm{\Lambda }}_x^2{Q}_{11}-9{\lambda }_h^2{\mathrm{\Lambda }}_x{Q}_{12}-12c_{\rho }/c_h^2{\omega }^2)Y=0& \hfill & \hfill \end{array}$

(15)

where λ _ρ = 2λ _h = −20, c _ρ = 7850, c _h = 0.01, Q₁₁ = Q₂₂ = Q₁₂/ν = Q₆₆/(1 − ν) = E/(1 − ν²), E = 200 GPa, and ν = 0.3 (x ∈ [0, 0.4] and y ∈ [0, 0.4]).

The selected exponential parameters in stiffness and density functions render the differential equation to have constant coefficients for which analytical solutions can be obtained. An analytical closed-form solution for the free-vibration frequencies can be obtained in this case and was used to benchmark the ASM.

The extended ASM’s stripwise exponential solution almost exactly matches the analytical solution for this exponentially graded plate, confirming the correctness of the method.

While the analytical solution seems to agree with ASM almost perfectly, the ANSYS solutions deviate significantly as the mode number increases. Table 3 presents the results obtained by MATLAB for ASM, ANSYS, and analytical method. An exponential change in plate density from 7850 kg/m³ to 1584.9 kg/m³ and thickness from 0.02 m to 0.009 m, without any stiffness property variation, resulted in a large discrepancy between ANSYS results and analytical results. The maximum percentage errors for ANSYS were 4.4%, 4.8%, and 4.0% for SSSS, SCSC, and SFSF cases, respectively. The same for ASM were 0.0%, 0.2%, and 0.0%.

Table 3.

Square plate with exponentially decreasing thickness and density: natural frequencies (Hz) for SSSS, SCSC, and SFSF. Analytical solution versus ASM (80 strips) and ANSYS (19440 elements).

BC	Method/Ref.	Mode shape
		(1,1)	(1,2)	(1,3)	(2,1)	(2,2)	(2,3)
SSSS	ANSYS	589.24	1479.0	2932.5	1476.0	2328.0	3749.5
	Analytical	606.86	1516.6	3023.8	1516.0	2418.1	3922.5
	ASM	606.86	1516.6	3023.8	1516.0	2418.1	3922.5
SCSC	ANSYS	889.74	2083.0	3801.2	1648.4	2905.9	4508.2
	Analytical	904.64	2138.2	3957.8	1685.8	2902.5	4735.6
	ASM	905.38	2141.1	3964.9	1686.7	2905.9	4743.8
SFSF	ANSYS	289.73	501.05	1103.6	1161.8	1410.2	2090.4
	Analytical	291.82	517.31	1142.4	1180.7	1445.3	2177.3
	ASM	291.82	517.31	1142.4	1180.7	1445.3	2177.3

Mode shapes related to the clamped edges in x direction are presented in Figure 5.

Figure 5.

Square plate with exponentially decreasing thickness and density: mode shapes for SCSC (ASM 80 strips).

3.3. Example 3

In this example, a square plate is considered with the same boundary conditions as in Example 2. However, in this case the bending rigidity and mass density vary with y in a power-law manner (a hypothetical material) chosen to produce a Cauchy–Euler equation upon separation. Specifically, the coefficients in the y-direction governing equation are defined as $Q_{11}=\frac{E}{1-{\nu }^2}{(c_L/y)}^6$ , $Q_{22}=\frac{E}{1-{\nu }^2}{(c_L/y)}^2$ , $Q_{12}=\nu \frac{E}{1-{\nu }^2}{(c_L/y)}^4$ , $Q_{66}=\frac{E}{1+\nu }{(c_L/y)}^4$ , while the density and thickness vary as ρ = c _ρ /y⁴ h = c _h y (x ∈ [0, 0.4] and y ∈ [1, 1.4]), respectively. The constant material properties are taken as E = 200 GPa, ν = 0.3, c _ρ = 7850 kg ⋅ m, c _h = 0.01, and c _L is a length parameter that ensures the dimensional integrity of the equation, in this example c _L = 1.

Under these spatially varying stiffness/mass/thickness assumptions, the equation of motion takes the form of a fourth-order Cauchy–Euler ODE, an equidimensional equation:

$\begin{array}{ccc}\hfill \frac{Y^{\prime \prime \prime \prime }}{y^2}+2\frac{Y^{\prime \prime \prime }}{y^3}-2c_L^2{\mathrm{\Lambda }}_x\frac{Y^{\prime \prime }}{y^4}+2c_L^2{\mathrm{\Lambda }}_x\frac{Y^{\prime }}{y^5}& \hfill & \hfill \\ \hfill +\left(c_L^4{\mathrm{\Lambda }}_x^2-2\nu c_L^2{\mathrm{\Lambda }}_x-12\frac{(1-{\nu }^2)c_{\rho }}{E{c}_L^2{c}_h^2}{\omega }^2\right)\frac{Y}{y^6}=0& \hfill & \hfill \end{array}$

(16)

which when multiplied by y⁶, yields the canonical form of the Cauchy–Euler equation. This type of ODE is solved by assuming a power-law solution Y(y) = y ^b , leading to an algebraic characteristic equation for b. The example provides a rigorous test: the ASM’s solution can be compared to the “exact” solution obtained.

In numerical experiments, the strip method achieved excellent agreement with the Cauchy–Euler benchmark solution. This demonstrates the extended ASM’s versatility in handling highly non-uniform plates—even when the spatial variation yields a non-constant coefficient ODE of the Cauchy–Euler type. Table 4 presents the results obtained by ASM and the analytical method. Mode shapes related to the clamped edges in x direction are presented in Figure 6.

Table 4.

Square plate with Cauchy–Euler coefficients in y: natural frequencies (Hz) for SSSS, SCSC, SFSF, and additional mixed sets. Analytical power-law solution versus ASM (80 strips).

BC	Method/Ref.	Mode shape
		(1,1)	(1,2)	(1,3)	(2,1)	(2,2)	(2,3)
SSSS	Analytical	362.88	998.64	2058.3	813.20	1448.8	2508.3
	ASM	362.88	998.62	2058.3	813.19	1448.7	2508.3
SCSC	Analytical	577.09	1440.4	2721.2	949.59	1814.3	3106.1
	ASM	577.39	1441.7	2724.5	949.97	1815.8	3109.8
SFSF	Analytical	145.83	276.57	707.18	590.27	754.56	1249.3
	ASM	145.83	276.57	707.18	590.28	754.56	1249.3
SSSC	Analytical	449.03	1198.5	2368.8	868.90	1612.3	2786.5
	ASM	449.22	1199.5	2371.5	869.17	1613.4	2789.5
SCSS	Analytical	459.50	1211.4	2382.5	873.79	1621.7	2798.0
	ASM	459.49	1211.3	2382.4	873.79	1621.7	2798.0
SSSF	Analytical	198.13	532.53	1263.8	650.34	1022.6	1771.6
	ASM	198.13	532.51	1263.7	650.34	1022.6	1771.6
SFSS	Analytical	179.60	517.70	1249.1	627.38	1005.7	1756.2
	ASM	179.60	517.71	1249.2	627.38	1005.7	1756.2
SCSF	Analytical	226.09	657.66	1502.9	665.24	1119.4	1975.2
	ASM	226.09	657.63	1502.8	665.24	1119.4	1975.1
SFSC	Analytical	201.42	631.41	1475.2	638.65	1095.2	1949.5
	ASM	201.45	631.76	1476.5	638.70	1095.8	1951.2

Figure 6.

Square plate with Cauchy–Euler coefficients in y: mode shapes for SCSC (ASM 80 strips).

4. Discussion of results

Example 1, which is analyzed in other studies, one using transfer matrix method (TMM), and one using Rayleigh-Ritz (RR), allows polynomial thickness variation. Results show that accurate frequencies are achieved with only about 9–10 ASM strips/elements (Table 2), whereas ANSYS requires a far denser (approximately 20000 elements) mesh to reach a comparable accuracy. The maximum discrepancy between ASM and ANSYS results were less than 0.081%. While the first and the second frequency of Rayleigh-Ritz example are within the results of other methods, the third frequency is slightly outside the range of other results. Whereas, ASM closely follows the FEM results. Even with two elements, where the flexural rigidity of the second strip (Y₂) is nearly twice of the flexural rigidity of the first strip (Y₁), the method yields satisfactory results.

Example 2 uses exponential thickness and density profiles so chosen as to yield a constant-coefficient ODE, for which an analytical solution is obtained, providing a benchmark case. ASM solutions were within a maximum error of 0.18% compared to the analytical solution.

Example 3 is constructed to produce a Cauchy–Euler-type equation, which again yielded an analytical solution for confirmation. Despite the strongly varying properties, ASM results were in excellent agreement with the analytical results (maximum error of 0.13%), showcasing the robustness of the method.

Compared to finite element analysis, the number of modes that can be obtained by ASM is not limited by the number of elements and the number of elements necessary for a certain accuracy is much less in ASM.

Furthermore, a great variety of spatially variable (1D) properties and thickness can be easily handled by the current method. Finally, ASM can also be utilized for free vibration of thin plates with loads (concentrated or distributed) or springs, as well as static solutions by setting ω = 0.

5. Conclusion

A semi-analytical extension of ASM applicable to thin plates with variable properties along one axis has been developed for free vibration analysis. Typical applications include tapered wing and rotor skins, bridge and ship panels with intentionally varied thickness, and functionally graded plate components designed to shift resonant frequencies, reduce weight, or redistribute stiffness in a prescribed manner. Plates with two opposite edges along the non-varying axis being simply-supported and the remaining boundary conditions being any combination of classical boundary conditions (clamped, simply-supported, free) are studied. In order to demonstrate the validity and accuracy of the method three numerical examples are presented.

• A semi-analytical extension of ASM was developed for free vibration analysis of thin plates with properties varying along one axis, under classical boundary conditions.

• The method was validated through three numerical examples (FEM, TMM, RR, and analytical) and showed excellent accuracy.

• In the polynomial-thickness case, ASM achieved accurate results with only about 9–10 elements, whereas ANSYS required a much denser mesh, demonstrating the computational efficiency of the method.

• The method remained accurate even for strongly varying thickness, density, and stiffness distributions, confirming its robustness.

• In addition to free vibration analysis, the proposed ASM formulation can be extended to plates with loads or springs, as well as to static analysis.