Vibration Analysis of Multilayered Quasicrystal Annular Plates, Cylindrical Shells, and Truncated Conical Shells Filled with Fluid

Abstract

An approach to estimate the dynamic characteristic of multilayered three-dimensional cubic quasicrystal cylindrical shells, annular plates, and truncated conical shells with different boundary conditions is presented. These investigated structures can be in a vacuum, totally filled with quiescent fluid, and subjected to internal flowing fluid where the fluid is incompressible and inviscid. The velocity potential, Bernoulli’s equation, and the impermeability condition have been applied to the shell–fluid interface to obtain an explicit expression, from which the fluid pressure can be converted into the coupled differential equations in terms of displacement functions. The state-space method is formulated to quasicrystal linear elastic theory to derive the state equations for the three structures along the radial direction. The mixed supported boundary conditions are represented by means of the differential quadrature technique and Fourier series expansions. A global propagator matrix, which connects the field variables at the internal interface to those at the external interface for the whole structure, is further completed by joint coupling matrices to overcome the numerical instabilities. Numerical examples show the correctness of the proposed method and the influence of the semi-vertical angle, different boundary conditions, and the fluid debit on the natural frequencies and mode shapes for various geometries and boundary conditions.

1. Introduction

Multilayered structures in solid mechanics couple the performances of two different materials by means of the stacking sequence, adhesive, and intermetallic surface deposition. Different from crystal and amorphous materials, the ordered and quasi-periodic atomic arrangement of quasicrystal (QC) materials enables them to exhibit some complex structures and special properties in theoretical analysis and experiments, especially high stiffness-to-weight ratios, high resistivity, and excellent thermal-electro-elastic behaviors [1,2], which can be used to improve the structure performances of multilayered plates and shells. For optimal design and fabrication of the multilayered QC structures with desired properties, it is necessary to thoroughly understand the mechanical behaviors of multi-field coupled QC materials in complex environments. QC atomized powders were sprayed on steel substrates to evaluate the mechanical properties of QC coatings [3]. Meanwhile, Ferreira et al. [4] investigated the wear resistance and low friction of functionally graded aluminum at high temperatures, where the QC approximant phase played a critical role. The results showed that the coated samples possessed high abrasion resistance, high micro-hardness, and excellent insulation. In addition, QC grains were developed into structural reinforcement phases and cast into alloys [5,6], from which the elongation, strength, and ductility of the specimens were increased. Due to the mechanical properties of multilayered composite structures being superior to those of components made of the individual material, they are useful in a tremendous number of different engineering applications such as aerospace, medical devices, and semi-conductors.

To prevent the initiation and propagation of defects in the multilayered plate/shell model design, the vibrational frequencies and stress distributions of the multilayered cylindrical shells, truncated conical shells, and annular plates need to be accurately predicted and evaluated [7,8,9]. Ten different shell theories were adopted to derive the semi-analysis solutions for cylindrical shells, and numerical examples indicated that the results based on Soedel and Kennard-Simplified theories were more accurate than those based on other theories [10]. Żur [11] presented the free vibration analysis of functionally graded circular plates with discrete elements by using the Neumann series method. Based on the discrete singular convolution method and differential quadrature method, Ersoy et al. [12] obtained the exact solutions for truncated conical shells, circular cylindrical shells, and annular plates. Eshaghi et al. [13] developed a numerical model of circular/annular plates under variable magnetic flux by using the Ritz method, from which the dynamic characteristics of the structure were obtained. Furthermore, they verified the effectiveness of this theoretical model by experimenting. The free vibration problems of determining the field variables which occurred in functionally graded annular plates with general boundary conditions [14,15] were derived, and the results showed that smaller gradient factors and elastic coefficients can retard the natural frequencies. Safarpour et al. [16] derived an exact solution of three-dimensional (3D) bending and frequency for functionally graded graphene-platelet-reinforced composite (FG-GPLRC) porous circular/annular plates with various boundary conditions. Meanwhile, the effects of geometric and material factors and different boundary conditions on the static and free vibration of FG-GPLRC were discussed using the theory of elasticity and differential quadrature method [17,18]. In addition, semi-analysis solutions for fiber metal annular plates with mixed boundary conditions [19] were developed to obtain the vibration frequencies. In the same way, truncated conical shells also have a wide range of engineering applications as the basic structure elements, and have attracted the attention of researchers. Kerboua and Lakis [20] obtained the free vibration analysis of discrete–continuous functionally graded circular plates via the Neumann series method. The first-order shell theory and the Sanders nonlinear kinematics equations [21] were proposed to improve the thermal and mechanical stability of functionally graded truncated conical shells.

Multi-physics phenomena between incompressible fluid flows and truncated conical shells may be encountered in practical engineering applications [22], such as oil transportation, urban pipe network laying, and fluid transmission. Zippo et al. [23] summarized and presented some results of an experimental campaign on the dynamic interactions between an elastic structure and a fluid. In theoretical analysis, Izyan et al. [24] presented the numerical solutions and numerical examples of frequency parameters for truncated conical shells filled with quiescent fluid using Love’s first approximation theory. Based on the Mindlin thick-shell theory, Hien et al. [25,26] derived the vibration responses of conical–cylindrical–conical shells containing fluid, to investigate the effect of fluid level, semi-vertex angles, and lamination sequences on three different joined cross-ply composite conical shells. Rahmanian et al. [27] obtained the natural frequencies of the coupled-field problem for multilayered truncated conical shells with different boundary conditions by using the same first-order shear deformation shell theory. In addition, for truncated conical shells with flowing liquid, Kerboua et al. [28] presented the explicit solution for fluid pressure according to Sanders’ thin-shell theory, and the numerical results were used to present how the geometry, stress stiffening effect, and hydrodynamic effect of the combined application of internal fluid and external flow influence the dynamic behavior of the combined shell. The nonlinear shell theory, along with Green’s strains, was proposed to analyze the dynamic instabilities of truncated conical shells subjected to flowing fluid [29], and some features showed that the appropriate dispersion of carbon nanotubes can improve the stability of the overall structure. In the literature mentioned above, the various approximation plate and shell theories were employed to analyze the dynamic responses of these structures. However, these approximate theories may have certain assumptions and be easily affected by the thickness of the structure, resulting in numerical instability. In addition, the traditional propagator matrix method in these studies cannot deal with the numerical instabilities in the case of large aspect ratios and high-order frequencies for QC plates/shells.

In this paper, the dynamic solutions for multilayered 3D cubic QC plates/shells subjected to flowing fluid are presented. The effects of different structure models, including the cylindrical shells, annular plates, and truncated conical shells, are taken into account. Governing equations are derived in the form of partial differential equations by means of the state-space method. According to the differential quadrature (DQ) techniques and Fourier series expressions, the extra modifications to these state equations are performed to satisfy the different boundary conditions. The velocity potential and Bernoulli’s equations are implemented to determine pressures from quiescent or flowing fluids acting on the truncated conical shells. Different from the traditional propagator matrix method, the new propagator relation is reestablished to resolve numerical instabilities in the case of large aspect ratios and high-order frequencies for QC plates/shells. One numerical example is directed to shed light on the effectiveness of the presented method. Other numerical examples are obtained to analyze the influence of the semi-vertical angle, different boundary conditions, and the fluid debit on the natural frequencies and mode shapes for the QC plates/shells containing the quiescent or flowing inviscid fluid. In addition, these numerical results can be emplaced as benchmarks for the design, numerical modeling, and simulation of 3D cubic QC plates/shells.

2. Theoretical Formulation

As shown in Figure 1, a multilayered 3D cubic QC truncated conical shell composed of M layers with perfect bonding is considered which has the length L = b − a, thickness h, and semi-vertical angle α. Each layer has a thickness of h_p = r_p − r_p−1, where r_p−1 and r_p denote the p-th layer of the inner and outer surfaces. It is clear that the inner and outer surfaces of this whole structure are r₀ and r_M, respectively. In a cylindrical coordinate system (r, θ, z), the mechanical behaviors of the cylindrical shell and annular plate can be obtained at the semi-vertical angles α = 0 and π/2, respectively.

2.1. Basic Equations

Based on the QC linear elasticity theory, the constitutive equations of 3D cubic QC materials [30] can be expressed as follows:

$$\begin{array}{l}{\sigma }_{rr}=C_{11}{\epsilon }_{rr}+C_{12}{\epsilon }_{\theta \theta }+C_{12}{\epsilon }_{zz}+R_1{w}_{rr}+R_2{w}_{\theta \theta }+R_2{w}_{zz},\\ {\sigma }_{\theta \theta }=C_{12}{\epsilon }_{rr}+C_{11}{\epsilon }_{\theta \theta }+C_{12}{\epsilon }_{zz}+R_2{w}_{rr}+R_1{w}_{\theta \theta }+R_2{w}_{zz},\\ {\sigma }_{zz}=C_{12}{\epsilon }_{rr}+C_{12}{\epsilon }_{\theta \theta }+C_{11}{\epsilon }_{zz}+R_2{w}_{rr}+R_2{w}_{\theta \theta }+R_1{w}_{zz},\\ {\sigma }_{\theta z}={\sigma }_{z\theta }=2C_{44}{\epsilon }_{z\theta }+2R_3{w}_{z\theta },\\ {\sigma }_{rz}={\sigma }_{zr}=2C_{44}{\epsilon }_{rz}+2R_3{w}_{rz},\\ {\sigma }_{r\theta }={\sigma }_{\theta r}=2C_{44}{\epsilon }_{r\theta }+2R_3{w}_{r\theta },\\ H_{rr}=R_1{\epsilon }_{rr}+R_2{\epsilon }_{\theta \theta }+R_2{\epsilon }_{zz}+K_{11}{w}_{rr}+K_{12}{w}_{\theta \theta }+K_{12}{w}_{zz},\\ H_{\theta \theta }=R_2{\epsilon }_{rr}+R_1{\epsilon }_{\theta \theta }+R_2{\epsilon }_{zz}+K_{12}{w}_{rr}+K_{11}{w}_{\theta \theta }+K_{12}{w}_{zz},\\ H_{zz}=R_2{\epsilon }_{rr}+R_2{\epsilon }_{\theta \theta }+R_1{\epsilon }_{zz}+K_{12}{w}_{rr}+K_{12}{w}_{\theta \theta }+K_{11}{w}_{zz},\\ H_{z\theta }=H_{\theta z}=2R_3{\epsilon }_{z\theta }+2K_{44}{w}_{z\theta },\\ H_{rz}=H_{zr}=2R_3{\epsilon }_{rz}+2K_{44}{w}_{rz},\\ H_{r\theta }=H_{\theta r}=2R_3{\epsilon }_{r\theta }+2K_{44}{w}_{r\theta },\end{array}$$

(1)

where σ_mn and H_mn (m, n = r, θ, z) represent the stress components in the phonon and phason fields, respectively. ε_mn and w_mn denote the phonon and phason strains, respectively. C₁₁, C₁₂, and C₄₄ are the phonon elastic constants, K₁₁, K₁₂, and K₄₄ denote the phason elastic constants, and R₁, R₂, and R₃ represent the phonon–phason elastic constants.

The relationship between the displacements and strains in phonon and phason fields [20,31] can be written as follows:

$$\begin{array}{l}{\epsilon }_{rr}={\displaystyle \frac{\partial u_r}{\partial r}}, {\epsilon }_{\theta \theta }=\tau \left({\displaystyle \frac{\partial u_{\theta }}{\partial \theta }}+u_r\cos \alpha +u_z\sin \alpha \right), {\epsilon }_{zz}={\displaystyle \frac{\partial u_z}{\partial z}},\\ {\epsilon }_{\theta z}={\epsilon }_{z\theta }={\displaystyle \frac{\tau }{2}}\left({\displaystyle \frac{\partial u_z}{\partial \theta }}-u_{\theta }\sin \alpha \right)+{\displaystyle \frac{1}{2}}{\displaystyle \frac{\partial u_{\theta }}{\partial z}}, {\epsilon }_{rz}={\epsilon }_{zr}={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\partial u_r}{\partial z}}+{\displaystyle \frac{\partial u_z}{\partial r}}\right), \\ {\epsilon }_{r\theta }={\epsilon }_{\theta r}={\displaystyle \frac{\tau }{2}}\left({\displaystyle \frac{\partial u_r}{\partial \theta }}-u_{\theta }\cos \alpha \right)+{\displaystyle \frac{1}{2}}{\displaystyle \frac{\partial u_{\theta }}{\partial r}},\\ w_{rr}={\displaystyle \frac{\partial w_r}{\partial r}}, w_{\theta \theta }=\tau \left({\displaystyle \frac{\partial w_{\theta }}{\partial \theta }}+w_r\cos \alpha +w_z\sin \alpha \right), w_{zz}={\displaystyle \frac{\partial w_z}{\partial z}},\\ w_{\theta z}=w_{z\theta }={\displaystyle \frac{\tau }{2}}\left({\displaystyle \frac{\partial w_z}{\partial \theta }}-w_{\theta }\sin \alpha \right)+{\displaystyle \frac{1}{2}}{\displaystyle \frac{\partial w_{\theta }}{\partial z}}, w_{rz}=w_{zr}={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\partial w_r}{\partial z}}+{\displaystyle \frac{\partial w_z}{\partial r}}\right), \\ w_{r\theta }=w_{\theta r}={\displaystyle \frac{\tau }{2}}\left({\displaystyle \frac{\partial w_r}{\partial \theta }}-w_{\theta }\cos \alpha \right)+{\displaystyle \frac{1}{2}}{\displaystyle \frac{\partial w_{\theta }}{\partial r}},\end{array}$$

(2)

where $\tau ={\displaystyle \frac{1}{r\cos \alpha +z\sin \alpha }}$ and u_m and w_m represent the phonon and phason displacements, respectively.

The body forces are ignored and the equations of motion for truncated conical shell [23,30,32] containing fluid can be defined as follows:

$$\begin{array}{l}{\displaystyle \frac{\partial {\sigma }_{rr}}{\partial r}}+\tau \left({\displaystyle \frac{\partial {\sigma }_{r\theta }}{\partial \theta }}+{\sigma }_{rz}\sin \alpha +{\sigma }_{rr}\cos \alpha -{\sigma }_{\theta \theta }\cos \alpha \right)+{\displaystyle \frac{\partial {\sigma }_{rz}}{\partial z}}=\rho {\displaystyle \frac{{\partial }^2{u}_r}{\partial t^2}}-p,\\ {\displaystyle \frac{\partial {\sigma }_{\theta r}}{\partial r}}+\tau \left({\displaystyle \frac{\partial {\sigma }_{\theta \theta }}{\partial \theta }}+2{\sigma }_{\theta z}\sin \alpha +2{\sigma }_{\theta r}\cos \alpha \right)+{\displaystyle \frac{\partial {\sigma }_{\theta z}}{\partial z}}=\rho {\displaystyle \frac{{\partial }^2{u}_{\theta }}{\partial t^2}},\\ {\displaystyle \frac{\partial {\sigma }_{zr}}{\partial r}}+\tau \left({\displaystyle \frac{\partial {\sigma }_{z\theta }}{\partial \theta }}\cos \alpha +{\sigma }_{zz}\sin \alpha -{\sigma }_{\theta \theta }\sin \alpha +{\sigma }_{zr}\cos \alpha \right)+{\displaystyle \frac{\partial {\sigma }_{zz}}{\partial z}}=\rho {\displaystyle \frac{{\partial }^2{u}_z}{\partial t^2}},\\ {\displaystyle \frac{\partial H_{rr}}{\partial r}}+\tau \left({\displaystyle \frac{\partial H_{r\theta }}{\partial \theta }}+H_{rz}\sin \alpha +H_{rr}\cos \alpha -H_{\theta \theta }\cos \alpha \right)+{\displaystyle \frac{\partial H_{rz}}{\partial z}}=\rho {\displaystyle \frac{{\partial }^2{w}_r}{\partial t^2}},\\ {\displaystyle \frac{\partial H_{\theta r}}{\partial r}}+\tau \left({\displaystyle \frac{\partial H_{\theta \theta }}{\partial \theta }}+2H_{\theta z}\sin \alpha +2H_{\theta r}\cos \alpha \right)+{\displaystyle \frac{\partial H_{\theta z}}{\partial z}}=\rho {\displaystyle \frac{{\partial }^2{w}_{\theta }}{\partial t^2}},\\ {\displaystyle \frac{\partial H_{zr}}{\partial r}}+\tau \left({\displaystyle \frac{\partial H_{z\theta }}{\partial \theta }}\cos \alpha +H_{zz}\sin \alpha -H_{\theta \theta }\sin \alpha +H_{zr}\cos \alpha \right)+{\displaystyle \frac{\partial H_{zz}}{\partial z}}=\rho {\displaystyle \frac{{\partial }^2{w}_z}{\partial t^2}},\end{array}$$

(3)

where p is the pressure exerted by the fluid on the shell wall; ρ denotes the mass density; t represents the time.

Based on a mixed formulation of solid mechanics, the state vector approach sets out-of-plane unknowns (u_r, u_θ, u_z, w_r, w_θ, w_z, σ_rr, σ_rθ, σ_rz, H_rr, H_rθ, H_rz) as basic variables, which are defined as primary variables A [33,34]. The state vector equations of the p-th layer can be established by combining the basic unknowns of the three governing equations in Equations (1)–(3), which are expressed as follows:

$${\displaystyle \frac{\partial A}{\partial r}}=MA=\left[\begin{array}{cc}{M}_1& M_2\\ M_3& M_4\end{array}\right]A,$$

(4)

where the submatrices M₁, M₂, M₃, and M₄ in Equation (4) can be found in Equations (A1)–(A4) of Appendix A, where some coefficients in Equations (A1)–(A4) and the following equations are listed in Equations (A5) and (A6).

The generalized displacement and stress solutions achieved by the Fourier trigonometric expansions [35,36] are proposed separately to solve the partial differential Equation (4) in the θ-direction as follows:

$$\left[\begin{array}{c}{u}_r\left(r,\theta ,z\right)\\ u_{\theta }\left(r,\theta ,z\right)\\ u_z\left(r,\theta ,z\right)\\ w_r\left(r,\theta ,z\right)\\ w_y\left(r,\theta ,z\right)\\ w_z\left(r,\theta ,z\right)\end{array}\right]={\displaystyle \sum _{n=1}^{\infty }\left[\begin{array}{c}{\tilde{u}}_r\left(r,z\right)\sin \left(n\theta \right)\\ {\tilde{u}}_{\theta }\left(r,z\right)\cos \left(n\theta \right)\\ {\tilde{u}}_z\left(r,z\right)\sin \left(n\theta \right)\\ {\tilde{w}}_r\left(r,z\right)\sin \left(n\theta \right)\\ {\tilde{w}}_{\theta }\left(r,z\right)\cos \left(n\theta \right)\\ {\tilde{w}}_z\left(r,z\right)\sin \left(n\theta \right)\end{array}\right]e^{\mathrm{i}\omega t}},$$

(5)

and

$$\left[\begin{array}{c}{\sigma }_{rr}\left(r,\theta ,z\right)\\ {\sigma }_{r\theta }\left(r,\theta ,z\right)\\ {\sigma }_{rz}\left(r,\theta ,z\right)\\ H_{rr}\left(r,\theta ,z\right)\\ H_{r\theta }\left(r,\theta ,z\right)\\ H_{rz}\left(r,\theta ,z\right)\end{array}\right]={\displaystyle \sum _{n=1}^{\infty }\left[\begin{array}{c}{\tilde{\sigma }}_{rr}\left(r,z\right)\sin \left(n\theta \right)\\ {\tilde{\sigma }}_{r\theta }\left(r,z\right)\cos \left(n\theta \right)\\ {\tilde{\sigma }}_{rz}\left(r,z\right)\sin \left(n\theta \right)\\ {\tilde{H}}_{rr}\left(r,z\right)\sin \left(n\theta \right)\\ {\tilde{H}}_{r\theta }\left(r,z\right)\cos \left(n\theta \right)\\ {\tilde{H}}_{rz}\left(r,z\right)\sin \left(n\theta \right)\end{array}\right]}{e}^{\mathrm{i}\omega t},$$

(6)

where n is the half wavenumber in the θ-direction; ω represents the angular frequency; The imaginary number i = $\sqrt{-1}$. It is difficult to deal with the partial differential equations of Equation (4) in the z-direction. Therefore, numerical manipulation, such as the DQ technique [37], is utilized to derive an exact solution for this problem. The mth-order derivative of an unknown function f(z) [14,38] at discrete point i can be written as follows:

$${\left.{\displaystyle \frac{{\partial }^{m}f\left(z_i\right)}{\partial z^m}}\right|}_{z=z_i}={\displaystyle \sum _{j=1}^N{g}_{ij}^{\left(m\right)}f\left(z_j\right)}\left(i=1,2,\dots ,N\right),$$

(7)

where $g_{ij}^{\left(m\right)}$ are weight coefficients.

The Chebyshev–Gauss–Lobatto grid space model [15,38] is employed to obtain the discrete points in the z-direction as follows:

$$z_i={\displaystyle \frac{b-a}{2}}\left[1-\cos \left({\displaystyle \frac{i-1}{N-1}}\pi \right)\right]+a.$$

(8)

2.2. Modeling of the Fluid–Shell Interaction

The dynamic response of the QC truncated conical shell may be influenced tremendously when subjected to flowing liquid. This liquid filling model is based on the following hypotheses [20,27]: (i) the internal fluid is inviscid and incompressible; (ii) the fluid velocity is constant; (iii) the vibration is linear. Based on these assumptions, the velocity potential ϕ satisfying the Laplace equation [39] is written as follows:

$${\nabla }^2\varphi ={\displaystyle \frac{{\partial }^2\varphi }{\partial z^2}}{\displaystyle \frac{1}{z^2{\sin }^2\xi }}{\displaystyle \frac{{\partial }^2\varphi }{\partial {\theta }^2}}+{\displaystyle \frac{2}{z}}{\displaystyle \frac{\partial \varphi }{\partial z}}+{\displaystyle \frac{1}{z^2}}{\displaystyle \frac{{\partial }^2\varphi }{\partial {\xi }^2}}+{\displaystyle \frac{1}{z^2\tan \xi }}{\displaystyle \frac{\partial \varphi }{\partial \xi }}=0,$$

(9)

where ξ (0 ≤ ξ ≤ α) is designated as the coordinate along the α-direction.

The relationship between the velocity potential and fluid velocity can be expressed as

$$V=\nabla \varphi ,$$

(10)

where V is the fluid velocity across a truncated conical shell section. In the absence of the net flow rate of fluid, the fluid velocity components V_θ, V_z, and V_α can be written as follows:

$$V_{\theta }={\displaystyle \frac{1}{z\sin \xi }}{\displaystyle \frac{\partial \varphi }{\partial \theta }}, V_z={\displaystyle \frac{Q}{z^2}}+{\displaystyle \frac{\partial \varphi }{\partial z}}, V_{\alpha }={\displaystyle \frac{1}{z}}{\displaystyle \frac{\partial \varphi }{\partial \xi }},$$

(11)

where Q is the fluid flowing with

$$Q=V_{ave}{z}^2=\mathrm{constant},$$

(12)

in which V_ave is the average fluid velocity.

The velocity potential ϕ [20] is designated as follows:

$$\varphi \left(\theta ,z,\xi ,t\right)=R\left(\xi \right)\phi \left(\theta ,z,t\right).$$

(13)

In order to satisfy the impermeability condition between the truncated conical shell and internal fluid, the fluid velocity on the internal surface must be consistent with the instantaneous rate of change of displacement as follows:

$${\left.{\left(V_{\alpha }\right)}_{\alpha =l}={\displaystyle \frac{1}{z}}{\displaystyle \frac{\partial \varphi }{\partial \xi }}\right|}_{\xi =l}={\left({\displaystyle \frac{\partial {\tilde{u}}_r}{\partial t}}+Q{\displaystyle \frac{\partial {\tilde{u}}_r}{\partial z}}\right)}_{\xi =l},$$

(14)

with l = α − tan(h/2z).

Substituting Equation (13) into Equation (14), the impermeability condition can be rewritten as

$$\phi \left(\theta ,z,t\right)={\displaystyle \frac{z}{R^{\prime }\left(l\right)}}{\left({\displaystyle \frac{\partial {\tilde{u}}_r}{\partial t}}+Q{\displaystyle \frac{\partial {\tilde{u}}_r}{\partial z}}\right)}_{\xi =l},$$

(15)

where the superscript “′” is the partial differentiation with respect to ξ. Then, combining Equation (13) with Equation (15), the representation of velocity potential ϕ is reestablished as follows:

$$\varphi \left(\theta ,z,\xi ,t\right)={\displaystyle \frac{zR\left(\xi \right)}{R^{\prime }\left(l\right)}}{\left({\displaystyle \frac{\partial {\tilde{u}}_r}{\partial t}}+Q{\displaystyle \frac{\partial {\tilde{u}}_r}{\partial z}}\right)}_{\xi =l}.$$

(16)

By plugging Equations (5), (6) and (16) into Equation (9), the function R(ξ) becomes

$$R^{″}\left(\xi \right)+{\displaystyle \frac{1}{\tan \xi }}{R}^{\prime }\left(\xi \right)-{\displaystyle \frac{n^2}{{\sin }^2\xi }}R\left(\xi \right)=0.$$

(17)

The solution of Equation (17) can be derived by using the Frobenius method [23] as follows:

$$R\left(\xi \right)=E{\xi }^n\left\{1+{\displaystyle \frac{n}{12}}{\xi }^2+{\displaystyle \frac{\left(5n+7\right)n}{1440}}{\xi }^4\right\},$$

(18)

where E is unknown.

Bernoulli’s equation [20] for the pressure of the liquid on the shell wall is given by

$$p=-{\rho }_f{\left({\displaystyle \frac{\partial \varphi }{\partial t}}+Q{\displaystyle \frac{\partial \varphi }{\partial z}}\right)}_{\xi =l},$$

(19)

where ρ_f is the fluid density.

Applying Equation (18) to Equation (13) and then combining it with Equation (19), the pressure p exerted by the fluid on the internal wall can be derived as follows:

$$p=-{\rho }_{f}k\left(Q^2{\displaystyle \frac{\partial {\tilde{u}}_r}{\partial z}}+Q^{2}z{\displaystyle \frac{{\partial }^2{\tilde{u}}_r}{\partial z^2}}+Q{\displaystyle \frac{\partial {\tilde{u}}_r}{\partial t}}+2Qz{\displaystyle \frac{{\partial }^2{\tilde{u}}_r}{\partial z\partial t}}+z{\displaystyle \frac{{\partial }^2{\tilde{u}}_r}{\partial t^2}}\right),$$

(20)

where k is a coefficient that is obtained from Equations (16) and (18). Applying the DQ method and Fourier trigonometric expansions to Equation (20), the pressure p can be rewritten as follows:

$$p_i=-{\rho }_{f}k\left(Q^2{\displaystyle \sum _{j=1}^N{X}_{ij}^{\left(1\right)}{\tilde{u}}_{rj}}+Q^2{z}_i{\displaystyle \sum _{j=1}^N{X}_{ij}^{\left(1\right)}{\tilde{u}}_{rj}}+Q\mathrm{i}\omega {\tilde{u}}_{ri}+2Q{z}_i\mathrm{i}\omega {\displaystyle \sum _{j=1}^N{X}_{ij}^{\left(1\right)}{\tilde{u}}_{rj}}-z_i{\omega }^2{\tilde{u}}_{ri}\right),$$

(21)

Thus, the first equation of motion in Equation (3) is related to phonon stresses and pressure.

2.3. Dispersion of State-Space Equations

The stresses and displacements in the θ-direction and z-direction are expressed in terms of Fourier trigonometric expansions and DQ regional discrete point expansions, respectively, from which the partial differential state equations in Equation (4) can be converted into ordinary differential equations as follows:

$${\displaystyle \frac{\mathrm{d}{\delta }^{\left(p\right)}}{\mathrm{d}r}}=T^{\left(p\right)}{\delta }^{\left(p\right)}={\left[\begin{array}{cc}{T}_1& T_2\\ T_3& T_4\end{array}\right]}^{(p)}{\delta }^{\left(p\right)},$$

(22)

where ${\delta }^{\left(p\right)}={\left[{\tilde{u}}_{ri},{\tilde{u}}_{\theta i},{\tilde{u}}_{zi},{\tilde{w}}_{ri},{\tilde{w}}_{\theta i},{\tilde{w}}_{zi},{\tilde{\sigma }}_{rri},{\tilde{\sigma }}_{r\theta i},{\tilde{\sigma }}_{rzi},{\tilde{H}}_{rri},{\tilde{H}}_{r\theta i},{\tilde{H}}_{rzi}\right]}^{\mathrm{T}}$; “T” represents the matrix or vector transpose; and the superscript “p” represents the state variables or coefficient matrix of the p-th layer. The submatrices T₁, T₂, T₃, and T₄ of the coefficient matrix T can be found in Equations (A7)–(A10) of Appendix A.

In this paper, simple supported (S) and clamped supported (C) boundary conditions are considered for multilayered cubic QC truncated conical shells. The mechanical boundary conditions at z = a and b can be written as follows [40]:

$$\mathrm{Simply} \mathrm{supported} (\mathrm{S}): {\tilde{u}}_{ri}={\tilde{u}}_{\theta i}={\tilde{w}}_{ri}={\tilde{w}}_{\theta i}=0, {\tilde{\sigma }}_{zzi}={\tilde{H}}_{zzi}=0,$$

(23)

$$\mathrm{Clamped} \mathrm{supported} \left(\mathrm{C}\right): {\tilde{u}}_{ri}={\tilde{u}}_{\theta i}={\tilde{u}}_{zi}={\tilde{w}}_{ri}={\tilde{w}}_{\theta i}={\tilde{w}}_{zi}=0.$$

(24)

For example, ‘CS’ denotes a shell with the clamped supported and simple supported boundary conditions at z = a and b, respectively. The semi-analysis solutions for QC truncated conical shells with four kinds of boundary conditions (SS, CS, SC, and CC) are obtained and the coefficient matrices of the state equations for shells with CS boundary conditions are listed in Equations (A13)–(A16) of Appendix A.

2.4. Vibration Analysis

Based on the theory of ordinary differential equations, the matrix exponential method can be utilized to derive the solution of Equation (22) [41] as follows:

$${\delta }^{\left(p\right)}\left(r\right)=\exp \left[T^{\left(p\right)}\left(r-r_{p-1}\right)\right]{\delta }^{\left(p\right)}\left(r_{p-1}\right)\left(r_{p-1}\le r\le r_p\right).$$

(25)

If r = r_p in Equation (25), this equation is in the following form:

$$\begin{gathered} {\delta }^{\left(p\right)}\left(r_p\right)=M^{\left(p\right)}{\delta }^{\left(p\right)}\left(r_{p-1}\right), \\ {M}^{\left(p\right)}=\exp \left[T^{\left(p\right)}\left(r-r_{p-1}\right)\right]. \end{gathered}$$

(26)

Equation (26) can be rewritten as follows:

$$\left[\begin{array}{c}{\delta }^{\left(p\right)+}\\ {\delta }^{\left(p\right)-}\end{array}\right]=\left[\begin{array}{c}{M}^{\left(p\right)}\\ I\end{array}\right]{\delta }^{\left(p\right)-}.$$

(27)

where the superscripts ‘+’ and ‘−’ are defined as variables on the upper surface and the lower surface of the p-th layer, respectively.

It is assumed that the field variables are continuous through the interface r_p as follows:

$${\delta }^{\left(p\right)+}={\delta }^{\left(p+1\right)-}.$$

(28)

According to the continuity interface conditions of r_p, Equation (28) can be redefined as follows:

$$J_c^{\left(p\right)}\left[\begin{array}{c}{\delta }^{\left(p+1\right)-}\\ {\delta }^{\left(p\right)+}\end{array}\right]=0\hspace{1em}\left(p=1,2,3,\dots ,M-1\right),$$

(29)

where $J_c^{\left(p\right)}$ = [I, −I] is the joint coupling matrix at the interface r_p.

Similarly, the interface relationship of the top and bottom surfaces for the QC truncated conical shell is

$$J_t{\delta }^{\left(M\right)+}=0, J_w{\delta }^{\left(1\right)-}=0,$$

(30)

where J_t = J_w = [I, 0].

Combining Equation (29) with (30), the interface condition of the integral shell is

$$J\delta =0,$$

(31)

where

$$\begin{array}{l}J=\mathrm{diag}\left[J_t, J_c, \cdots , J_c, J_w\right], \\ \delta ={\left[{\left({\delta }^{\left(M\right)+}\right)}^{\mathrm{T}},\dots ,{\left(\begin{array}{c}{\delta }^{\left(p+1\right)-}\\ {\delta }^{\left(p\right)+}\end{array}\right)}^{\mathrm{T}},\dots , {\left({\delta }^{\left(1\right)-}\right)}^{\mathrm{T}}\right]}^{\mathrm{T}}.\end{array}$$

(32)

The whole structure can be defined in the same form as Equation (31) as follows:

$$\delta =\Pi {\delta }^{-},$$

(33)

where

$$\begin{array}{l}\Pi =\mathrm{diag}\left[\left(\begin{array}{c}{M}^{\left(M\right)}\\ I\end{array}\right), \left(\begin{array}{c}{M}^{\left(M-1\right)}\\ I\end{array}\right), \cdots , \left(\begin{array}{c}{M}^{\left(2\right)}\\ I\end{array}\right), \left(\begin{array}{c}{M}^{\left(1\right)}\\ I\end{array}\right)\right],\\ {\delta }^{-}={\left[\left({\delta }^{\left(M\right)-}\right), \left({\delta }^{\left(M-1\right)-}\right), \cdots , \left({\delta }^{\left(2\right)-}\right), \left({\delta }^{\left(1\right)-}\right)\right]}^{\mathrm{T}}.\end{array}$$

(34)

Substituting Equation (33) into Equation (31) results in

$$J\Pi {\delta }^{-}=0.$$

(35)

Equation (35) can be rewritten as follows:

$$\left|J\Pi \right|=0.$$

(36)

By solving Equation (36), the vibration solutions for the whole multilayered QC structure are obtained. As a summary, Figure 2 shows the solution procedure utilized in this study.

3. Numerical Results and Discussion

Illustrative examples of the state-space-based differential quadrature method (SS-DQM) for the multi-physics problems between incompressible fluid flows and truncated conical shells are presented. The aspect ratios for these structures are set as a/b = 0.4 and h/b = 0.03. The cylindrical shells, annular plates, and truncated conical shells are composed of two 3D cubic QC materials, and they are abbreviated as QC1 and QC2. Meanwhile, the material properties [1,30] of the constituents are listed in

Table 1. Material properties (C, K, and R in 10⁹ N/m² and ρ in Kg/m³).

	C11	C12	C44	K11	K12	K44	R1	R2	R3	ρ
QC1	234.33	57.41	70.19	122	24	12	8.846	4.578	3.125	4505
QC2	112.1	60.3	32.8	60	20	10	5	−2	7	5300

, and these material coefficients completely match the elastic deformation energy density of the QC materials and satisfy the positive definiteness of the QC materials [42,43]. The structure (called QC1/QC2/QC1), composed of three single layers of the same thickness, is taken in the following text. The half wavenumber in Equations (5) and (6) is taken as n = 1, and the discrete point is set to N = 13. In addition, the dimensionless quantities [44] are employed,

$${\sigma }_{mn}^{*}={\displaystyle \frac{{\sigma }_{mn}}{C_{\max }}}, H_{mn}^{*}={\displaystyle \frac{H_{mn}}{C_{\max }}}, u_m^{*}={\displaystyle \frac{u_m}{b}}, w_m^{*}={\displaystyle \frac{w_m}{b}}, {\tau }^{*}={\displaystyle \frac{\tau }{b}}, z^{*}={\displaystyle \frac{z}{b}}, r^{*}={\displaystyle \frac{r}{b}}, \mathsf{\Omega }=\omega b\sqrt{{\rho }_{\max }/C_{\max }},$$

(37)

where C_max and ρ_max denote the maximum values of the phonon elastic constant and mass density for QC multilayered structures, respectively.

3.1. Convergence and Validation Studies

To prove the formula and program in this paper, the numerical results of the QC model achieved by the proposed method are compared with those achieved by the state-space method (SSM) [33]. The first seven order dimensionless natural frequencies Ω of the QC cylindrical shells (α = 0) with boundary conditions SS are presented in

Table 2. Comparisons of SS-DQM results for the first seven non-dimensional natural frequencies Ω with the SSM results for QC cylindrical shells.

	Mode 1	Mode 2	Mode 3	Mode 4	Mode 5	Mode 6	Mode 7
Result by SSM	0.50665	0.75775	6.26025	12.88351	18.83285	23.90636	26.45587
Present result	0.50665	0.75775	6.26026	12.88351	18.83285	23.90637	26.45588

. It is noted that the first five order modes are the flexural frequencies, but Modes 6 and 7 are longitudinal frequencies. The maximal relative error of the frequency values Ω obtained from the present and SSM is less than 0.01%.

Similarly, Figure 3a shows that the variations of the first-order normalized mode shapes of the displacements and stresses along the r-direction are obtained using the present method and SSM. The distributions of u_r and w_r are flexural mode shapes. The results obtained by the two methods are consistent. This feature indicates that the method in this paper has high precision and a good convergence and the solutions are stable. As mentioned in the introduction, research works for the vibration analysis of empty and fluid-filled QC conical shells are encountered less frequently. Because cylindrical shells can be viewed as a special case of conical vessels with equal top and bottom radii, both Ω and mode shapes are found to have good consistency by using the present method and SSM, which can be used as an example to verify the correctness of formulas and procedures for QC conical shells.

To ensure that the present method is valid and accurate, a 3D cubic QC sandwich plate with simply supported boundary conditions is considered. The material coefficients, the shape and size, and the expression form of the field variables are consistent with those in Ref. [30]. In addition, the interface coefficient is taken as δ = 0. The numerical results for the QC1/QC2/QC1 plate are calculated by using the SS-DQM to compare with results of Ref. [30] in Figure 6a–d. The first five order dimensionless natural frequencies Ω of the QC plate are presented in

Table 3. Comparisons of SS-DQM results for first five non-dimensional natural frequencies Ω with the exact solutions of the pseudo-Stroh formalism of Ref. [30].

	Mode 1	Mode 2	Mode 3	Mode 4	Mode 5
Result by Ref. [30]	1.47999	2.38416	4.27451	8.32909	9.52440
Present result	1.47999	2.38417	4.27453	8.32911	9.52443

. Meanwhile, Figure 3b presents the variation of the phonon and phason displacements u_x, w_x, u_z, and w_z. The solutions obtained for Ω and normalized displacement mode shapes by using the SS-DQM are in good agreement with the results of Ref. [30] achieved by using pseudo-Stroh formalism. Thus, the method in this paper has high precision and good convergence.

3.2. QC Annular Plates with Different Boundary Conditions

In this part, the solutions for the QC annular plates (α = π/2) with the boundary conditions SS, CS, SC, and CC are presented. The first nine order dimensionless natural frequencies Ω of the QC annular plates without fluid are presented in

Table 4. The first nine non-dimensional natural frequencies Ω of QC annular plates with different boundary conditions.

	Mode 1	Mode 2	Mode 3	Mode 4	Mode 5	Mode 6	Mode 7	Mode 8	Mode 9
SS	0.18006	0.25376	0.66903	0.71981	0.95204	1.16316	1.40024	1.69379	2.05038
CS	0.25230	0.35814	0.79363	1.15242	1.30192	1.54856	1.96985	2.17479	2.32465
SC	0.28126	0.39808	0.82056	1.19096	1.52124	1.57224	2.23517	2.36010	2.53888
CC	0.36982	0.53056	0.94874	1.40814	1.54224	1.71593	2.57940	2.63754	2.78044

to show the effect of the simple supported and clamped supported boundary conditions on Ω. Every order Ω increases in the order of SS, CS, SC, and CC for a fixed value of h/L. Different from QC cylindrical shells with the boundary condition SS, Modes 1, 2, 3, 5, 7, and 9 of the annular plates are flexural frequencies, but Modes 4, 6, and 8 are longitudinal frequencies.

The radial coordinates of data points in the resulting graph of this section are fixed at z/b = 0.7. The ring coordinates of data points are all basic solutions, and the angle variation can be multiplied by the sine-cosine factor from Equations (5) and (6). Figure 4 shows the variations of the first-order mode shapes of the phonon and phason displacements for QC annular plates with four different boundary conditions along the r-direction. The distributions of $u_r^{*}$ and $w_r^{*}$ (Figure 4a,c) are symmetrical, but the distributions of $u_{\theta }^{*}$ and $w_z^{*}$ (Figure 4b,d) are anti-symmetrical at r/b = 0.045, so these curves are flexural mode shapes. Moreover, the values of the mode shapes for $u_r^{*}$ and $u_{\theta }^{*}$ at r/b = 0.03 and 0.06 decrease in the order of SC, CC, SS, and CS for a fixed value of h/L. In addition, it is observed that the phason mode shapes $w_r^{*}$ and $w_z^{*}$ are sensitive to the different boundary conditions.

Figure 5 shows the variations of the first-order mode shapes of the phonon and phason stresses for QC annular plates with four different boundary conditions along the r-direction. The ${\mathsf{\sigma }}_{\mathit{rr}}^{*}$ (Figure 5a) of QC annular plates with boundary conditions CS and CC experience a large slope jump at the two interfaces based on the QC linear elastic theory. In addition, increasing the phonon material properties of QC2 helps reduce this jump. The ${\mathsf{\sigma }}_{r\theta }^{*}$ (Figure 5b) of the plates with boundary conditions CS exhibits a different response from plates with other boundary conditions. The maximal values of $H_{\mathit{rr}}^{*}$ and $H_{\mathit{rz}}^{*}$ (Figure 5c,d) at r/b = 0.03 and 0.06 increase in the order of SS, SC, CS, and CC for a fixed value of h/L. These features show that the clamped boundary condition increases the stiffness of the plate, resulting in changes in the frequencies and mode shapes.

3.3. Dynamic Behavior of the Empty and Fluid-Filled Shell

In this part, the solutions for the QC truncated conical shells with quiescent fluid are presented. The first five order dimensionless natural frequencies Ω of the QC shells with boundary condition SS are listed in

Table 5. The first five non-dimensional natural frequencies Ω of QC truncated conical shells with different semi-vertical angles.

	Mode 1	Mode 2	Mode 3	Mode 4	Mode 5
α0 = π/6	0.55916	0.79145	1.13472	1.28815	1.60288
α = π/6	0.55959	0.79158	1.15277	1.28852	1.63267
α = π/4	0.38664	0.71957	1.01077	1.41469	1.60849
α = π/3	0.28030	0.47297	0.82724	1.21100	1.36207

, in which α₀ is defined as the shell without fluid. The natural frequencies Ω of shells containing quiescent fluid are greater than that of the shell without fluid at the same α. This feature shows that quiescent fluid improves the overall stiffness of these shells. Meanwhile, as α increases, Ω decreases, except for in Mode 4, indicating that the device composed of QC truncated conical shells would be more susceptible to external motivation with an increasing α. In other words, the QC structures with a large α are more likely to cause structural resonance under external motivation.

The contour plots of the dimensionless phonon displacement mode shapes $u_r^{*}$ and $u_{\theta }^{*}$ for QC truncated conical shells on the r-z plane are presented in Figure 6 and Figure 7. The distributions of $u_r^{*}$ (Figure 6a–d) for the QC shell are all asymmetrical at z/b = 0.7 mm and almost constant along the r-direction. The maximum value of $u_r^{*}$ becomes greater with the increase in α. The results in Figure 6 and Figure 7 confirm that the phonon displacement boundary conditions in Equation (27) are precisely satisfied, partially verifying the accuracy of the present methods. The distributions of $u_{\theta }^{*}$ (Figure 7a–d) are nonlinear, from which the first-order, second-order, and higher-order plate and shell theories are imprecise.

The r-z plane distributions of the dimensionless phason displacement mode shapes $w_{\theta }^{*}$ and $w_z^{*}$ for QC truncated conical shells are plotted in Figure 8 and Figure 9. $w_{\theta }^{*}$ (Figure 8a–d) experiences a large slope jump in the r-z plane for the classical elasticity case, and increasing α changes the direction of the jump. It is seen that the distribution of $w_z^{*}$ in Figure 9a–d is a polygonal line while that of $w_{\theta }^{*}$ is not, and $w_{\theta }^{*}$ and $w_z^{*}$ show different behaviors in the r-z plane. This difference is mainly due to differences in material properties since the two in-plane mode shapes of the QC shell should have the same distribution curves in a physical sense. $w_{\theta }^{*}$ and $w_z^{*}$ remain continuous at the interfaces at r/b = 0.04 and 0.05. The distributions of $u_{\theta }^{*}$ and $w_{\theta }^{*}$ are similar in the r-z plane. This feature indicates that the phonon and phason elementary excitations in QCs both need to be analyzed.

3.4. QC Truncated Conical Shells Containing Flowing Fluid

In this part, the solutions for the interaction between incompressible fluid flow and QC truncated conical shells with boundary conditions CS are presented. The centrifugal, Coriolis, and inertial forces will affect the stability of the whole structure when the incompressible fluid is flowing through the shell. The material and physical properties of the whole model are α = π/12 and ρ_f = 1000 kg/m³. The first five order dimensionless natural frequencies Ω of the QC shells with different fluid debits Q are listed in

Table 6. The first nine non-dimensional natural frequencies Ω of QC truncated conical shells subjected to flowing fluids.

	Mode 1	Mode 2	Mode 3	Mode 4	Mode 5	Mode 6	Mode 7	Mode 8	Mode 9
Q = 1 m3/s	0.83832	1.56236	1.75879	2.35217	2.42893	3.26236	3.84708	4.03837	4.17203
Q = 2 m3/s	0.83946	1.56357	1.75968	2.41162	2.43311	3.26878	3.85937	4.18399	4.58554
Q = 3 m3/s	0.83971	1.56387	1.75992	2.42240	2.43481	3.27066	3.85993	4.18815	4.59383
Q = 4 m3/s	0.83980	1.56399	1.76003	2.42608	2.43588	3.27141	3.86014	4.18988	4.59572
Q = 5 m3/s	0.83985	1.56405	1.76010	2.42784	2.43669	3.27177	3.86026	4.19075	4.59665

. These frequencies Ω are increased by increasing the mean fluid velocity due to Q being a constant related to the flow rate. Meanwhile, these Ω values show that the model composed of fluid flow and QC shells can be accurately calculated at a small α.

Figure 10 and Figure 11 show the variation of the first-order dimensionless displacement and stress mode shapes of QC truncated conical shells with different Q values, respectively. An increase in Q results in a decrease in the magnitude of $u_r^{*}$ and $u_{\theta }^{*}$ (Figure 10a,b) for QC shells. The decreasing of the values of the mode shapes is fast in the range of 1 m³/s ≤ Q ≤ 2 m³/s. It is noted that $u_r^{*}$ and $u_{\theta }^{*}$ are continuous across the interfaces for the classical elasticity case. However, $u_z^{*}$ (Figure 10c) shows different behavior, especially in the middle layers of the QC shell. $u_z^{*}$ experiences a slope jump at r/b = 0.04 and 0.05. In addition, Q has little effect on $w_r^{*}$ (Figure 10d), including $w_{\theta }^{*}$, and $w_z^{*}$. The maximum values of ${\mathsf{\sigma }}_{\mathit{rr}}^{*}$ and ${\mathsf{\sigma }}_{\mathit{rz}}^{*}$ (Figure 11a,c) occur at the interfaces between layers. Different from the distribution of $u_{\theta }^{*}$, ${\mathsf{\sigma }}_{r\theta }^{*}$ (Figure 11b) is insensitive to Q. Similarly, $H_{\mathit{rz}}^{*}$ (Figure 11d) is as insensitive to Q as the phason displacement mode shapes. This feature shows the weak coupling between phason and velocity potential.

4. Conclusions

This paper presents a semi-analysis method to evaluate the natural frequencies and mode shapes of 3D cubic QC cylindrical shells, annular plates, and truncated conical shells under vacuum and fluid. Since no assumptions on stresses and displacements have been employed, the present framework can be applied to an arbitrary thickness of plate/shell structures. Some significant features are listed below.

• The global propagator relation is reestablished to solve numerical instabilities in the case of large discrete point numbers and high frequencies for three QC structures. Some cases such as one-dimensional and two-dimensional QC plates/shells could be investigated according to the present methods.
• Quiescent fluid improves the overall stiffness of the truncated conical shell. This shell with the higher semi-vertical angle α leads to smaller natural frequencies and phonon displacement mode shapes and larger phason displacement mode shapes. Meanwhile, QC truncated conical shells are more susceptible to external excitation with the increase in semi-vertical angle α.
• Every order of natural frequencies increases in the order of SS, CS, SC, and CC for a fixed value of h/L. Higher fluid debits Q lead to increased frequencies for the shells across all boundary conditions. All phason variables are insensitive to fluid debits Q.

The methods and numerical results in this paper can be utilized to validate the accuracy of other numerical methods and serve for the analysis and design of intelligent QC material plates/shells.